diff --git a/test/Jacobian/FunctionCalls.C b/test/Jacobian/FunctionCalls.C
index da7d5420c..bd6f187f1 100644
--- a/test/Jacobian/FunctionCalls.C
+++ b/test/Jacobian/FunctionCalls.C
@@ -6,38 +6,25 @@
 #include <cmath>
 #include "clad/Differentiator/Differentiator.h"
 
-double outputs[4], results[4];
+double outputs[4];
+clad::matrix<double> results(2, 2);
 
 void fn1(double i, double j, double* output) {
   output[0] = std::pow(i, j);
   output[1] = std::pow(j, i);
 }
 
-// CHECK: void fn1_jac(double i, double j, double *output, double *jacobianMatrix) {
-// CHECK-NEXT:     double _t0 = output[0];
-// CHECK-NEXT:     output[0] = std::pow(i, j);
-// CHECK-NEXT:     double _t1 = output[1];
-// CHECK-NEXT:     output[1] = std::pow(j, i);
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             double _r2 = 0.;
-// CHECK-NEXT:             double _r3 = 0.;
-// CHECK-NEXT:             clad::custom_derivatives::pow_pullback(j, i, 1, &_r2, &_r3);
-// CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += _r2;
-// CHECK-NEXT:             jacobianMatrix[{{2U|2UL|2ULL}}] += _r3;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         output[1] = _t1;
-// CHECK-NEXT:     }
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             double _r0 = 0.;
-// CHECK-NEXT:             double _r1 = 0.;
-// CHECK-NEXT:             clad::custom_derivatives::pow_pullback(i, j, 1, &_r0, &_r1);
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += _r0;
-// CHECK-NEXT:             jacobianMatrix[{{1U|1UL|1ULL}}] += _r1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         output[0] = _t0;
-// CHECK-NEXT:     }
+// CHECK: void fn1_jac(double i, double j, double *output, clad::array_ref<clad::matrix<double> > _d_vector_output) {
+// CHECK-NEXT:     unsigned long indepVarCount = * _d_vector_output.rows() + {{2U|2UL|2ULL}};
+// CHECK-NEXT:     clad::array<double> _d_vector_i = clad::one_hot_vector(indepVarCount, {{0U|0UL|0ULL}});
+// CHECK-NEXT:     clad::array<double> _d_vector_j = clad::one_hot_vector(indepVarCount, {{1U|1UL|1ULL}});
+// CHECK-NEXT:     * _d_vector_output = clad::identity_matrix(* _d_vector_output.rows(), indepVarCount, {{2U|2UL|2ULL}});
+// CHECK-NEXT:     {{.*}} _t0 = clad::custom_derivatives::pow_pushforward(i, j, _d_vector_i, _d_vector_j);
+// CHECK-NEXT:     * _d_vector_output[0] = _t0.pushforward;
+// CHECK-NEXT:     output[0] = _t0.value;
+// CHECK-NEXT:     {{.*}} _t1 = clad::custom_derivatives::pow_pushforward(j, i, _d_vector_j, _d_vector_i);
+// CHECK-NEXT:     * _d_vector_output[1] = _t1.pushforward;
+// CHECK-NEXT:     output[1] = _t1.value;
 // CHECK-NEXT: }
 
 #define INIT(F) auto d_##F = clad::jacobian(F);
@@ -47,17 +34,16 @@ void fn1(double i, double j, double* output) {
 template <unsigned numOfOutputs, typename Fn, typename... Args>
 void test(Fn derivedFn, Args... args) {
   unsigned numOfParameters = sizeof...(args);
-  unsigned numOfResults = numOfOutputs * numOfParameters;
   for (unsigned i = 0; i < numOfOutputs; ++i)
     outputs[i] = 0;
-  for (unsigned i = 0; i < numOfResults; ++i)
-    results[i] = 0;
-  derivedFn.execute(args..., outputs, results);
+  derivedFn.execute(args..., outputs, &results);
   printf("{");
-  for (unsigned i = 0; i < numOfResults; ++i) {
-    printf("%.2f", results[i]);
-    if (i != numOfResults - 1)
-      printf(", ");
+  for (unsigned i = 0; i < numOfOutputs; ++i) {
+    for (unsigned j = 0; j < numOfParameters; ++j) {
+      printf("%.2f", results[i][j]);
+      if (i != numOfOutputs - 1 || j != numOfParameters - 1)
+        printf(", ");
+    }
   }
   printf("}\n");
 }
diff --git a/test/Jacobian/Functors.C b/test/Jacobian/Functors.C
index cf04e1291..f48599364 100644
--- a/test/Jacobian/Functors.C
+++ b/test/Jacobian/Functors.C
@@ -16,27 +16,21 @@ struct Experiment {
     x = val;
   }
 
-  // CHECK: void operator_call_jac(double i, double j, double *output, double *jacobianMatrix) {
-  // CHECK-NEXT:     double _t0 = output[0];
-  // CHECK-NEXT:     output[0] = this->x * i * i * j;
-  // CHECK-NEXT:     double _t1 = output[1];
-  // CHECK-NEXT:     output[1] = this->y * i * j * j;
-  // CHECK-NEXT:     {
-  // CHECK-NEXT:         {
-  // CHECK-NEXT:             jacobianMatrix[{{2U|2UL|2ULL}}] += this->y * 1 * j * j;
-  // CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += this->y * i * 1 * j;
-  // CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += this->y * i * j * 1;
-  // CHECK-NEXT:         }
-  // CHECK-NEXT:         output[1] = _t1;
-  // CHECK-NEXT:     }
-  // CHECK-NEXT:     {
-  // CHECK-NEXT:         {
-  // CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += this->x * 1 * j * i;
-  // CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += this->x * i * 1 * j;
-  // CHECK-NEXT:             jacobianMatrix[{{1U|1UL|1ULL}}] += this->x * i * i * 1;
-  // CHECK-NEXT:         }
-  // CHECK-NEXT:         output[0] = _t0;
-  // CHECK-NEXT:     }
+  // CHECK: void operator_call_jac(double i, double j, double *output, clad::array_ref<clad::matrix<double> > _d_vector_output) {
+  // CHECK-NEXT:     unsigned long indepVarCount = * _d_vector_output.rows() + {{2U|2UL|2ULL}};
+  // CHECK-NEXT:     clad::array<double> _d_vector_i = clad::one_hot_vector(indepVarCount, {{0U|0UL|0ULL}});
+  // CHECK-NEXT:     clad::array<double> _d_vector_j = clad::one_hot_vector(indepVarCount, {{1U|1UL|1ULL}});
+  // CHECK-NEXT:     * _d_vector_output = clad::identity_matrix(* _d_vector_output.rows(), indepVarCount, {{2U|2UL|2ULL}});
+  // CHECK-NEXT:     double &_t0 = this->x;
+  // CHECK-NEXT:     double _t1 = _t0 * i;
+  // CHECK-NEXT:     double _t2 = _t1 * i;
+  // CHECK-NEXT:     * _d_vector_output[0] = ((0 * i + _t0 * _d_vector_i) * i + _t1 * _d_vector_i) * j + _t2 * _d_vector_j;
+  // CHECK-NEXT:     output[0] = _t2 * j;
+  // CHECK-NEXT:     double &_t3 = this->y;
+  // CHECK-NEXT:     double _t4 = _t3 * i;
+  // CHECK-NEXT:     double _t5 = _t4 * j;
+  // CHECK-NEXT:     * _d_vector_output[1] = ((0 * i + _t3 * _d_vector_i) * j + _t4 * _d_vector_j) * j + _t5 * _d_vector_j;
+  // CHECK-NEXT:     output[1] = _t5 * j;
   // CHECK-NEXT: }
 
 };
@@ -52,27 +46,21 @@ struct ExperimentConst {
     x = val;
   }
 
-  // CHECK: void operator_call_jac(double i, double j, double *output, double *jacobianMatrix) const {
-  // CHECK-NEXT:     double _t0 = output[0];
-  // CHECK-NEXT:     output[0] = this->x * i * i * j;
-  // CHECK-NEXT:     double _t1 = output[1];
-  // CHECK-NEXT:     output[1] = this->y * i * j * j;
-  // CHECK-NEXT:     {
-  // CHECK-NEXT:         {
-  // CHECK-NEXT:             jacobianMatrix[{{2U|2UL|2ULL}}] += this->y * 1 * j * j;
-  // CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += this->y * i * 1 * j;
-  // CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += this->y * i * j * 1;
-  // CHECK-NEXT:         }
-  // CHECK-NEXT:         output[1] = _t1;
-  // CHECK-NEXT:     }
-  // CHECK-NEXT:     {
-  // CHECK-NEXT:         {
-  // CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += this->x * 1 * j * i;
-  // CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += this->x * i * 1 * j;
-  // CHECK-NEXT:             jacobianMatrix[{{1U|1UL|1ULL}}] += this->x * i * i * 1;
-  // CHECK-NEXT:         }
-  // CHECK-NEXT:         output[0] = _t0;
-  // CHECK-NEXT:     }
+  // CHECK: void operator_call_jac(double i, double j, double *output, clad::array_ref<clad::matrix<double> > _d_vector_output) const {
+  // CHECK-NEXT:     unsigned long indepVarCount = * _d_vector_output.rows() + {{2U|2UL|2ULL}};
+  // CHECK-NEXT:     clad::array<double> _d_vector_i = clad::one_hot_vector(indepVarCount, {{0U|0UL|0ULL}});
+  // CHECK-NEXT:     clad::array<double> _d_vector_j = clad::one_hot_vector(indepVarCount, {{1U|1UL|1ULL}});
+  // CHECK-NEXT:     * _d_vector_output = clad::identity_matrix(* _d_vector_output.rows(), indepVarCount, {{2U|2UL|2ULL}});
+  // CHECK-NEXT:     double &_t0 = this->x;
+  // CHECK-NEXT:     double _t1 = _t0 * i;
+  // CHECK-NEXT:     double _t2 = _t1 * i;
+  // CHECK-NEXT:     * _d_vector_output[0] = ((0 * i + _t0 * _d_vector_i) * i + _t1 * _d_vector_i) * j + _t2 * _d_vector_j;
+  // CHECK-NEXT:     output[0] = _t2 * j;
+  // CHECK-NEXT:     double &_t3 = this->y;
+  // CHECK-NEXT:     double _t4 = _t3 * i;
+  // CHECK-NEXT:     double _t5 = _t4 * j;
+  // CHECK-NEXT:     * _d_vector_output[1] = ((0 * i + _t3 * _d_vector_i) * j + _t4 * _d_vector_j) * j + _t5 * _d_vector_j;
+  // CHECK-NEXT:     output[1] = _t5 * j;
   // CHECK-NEXT: }
 
 };
@@ -88,29 +76,21 @@ struct ExperimentVolatile {
     x = val;
   }
 
-  // CHECK: void operator_call_jac(double i, double j, double *output, double *jacobianMatrix) volatile {
-  // CHECK-NEXT:     double _t0 = this->x * i;
-  // CHECK-NEXT:     double _t1 = output[0];
-  // CHECK-NEXT:     output[0] = this->x * i * i * j;
-  // CHECK-NEXT:     double _t2 = this->y * i;
-  // CHECK-NEXT:     double _t3 = output[1];
-  // CHECK-NEXT:     output[1] = this->y * i * j * j;
-  // CHECK-NEXT:     {
-  // CHECK-NEXT:         {
-  // CHECK-NEXT:             jacobianMatrix[{{2U|2UL|2ULL}}] += this->y * 1 * j * j;
-  // CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += _t2 * 1 * j;
-  // CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += _t2 * j * 1;
-  // CHECK-NEXT:         }
-  // CHECK-NEXT:         output[1] = _t3;
-  // CHECK-NEXT:     }
-  // CHECK-NEXT:     {
-  // CHECK-NEXT:         {
-  // CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += this->x * 1 * j * i;
-  // CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += _t0 * 1 * j;
-  // CHECK-NEXT:             jacobianMatrix[{{1U|1UL|1ULL}}] += _t0 * i * 1;
-  // CHECK-NEXT:         }
-  // CHECK-NEXT:         output[0] = _t1;
-  // CHECK-NEXT:     }
+  // CHECK: void operator_call_jac(double i, double j, double *output, clad::array_ref<clad::matrix<double> > _d_vector_output) volatile {
+  // CHECK-NEXT:     unsigned long indepVarCount = * _d_vector_output.rows() + {{2U|2UL|2ULL}};
+  // CHECK-NEXT:     clad::array<double> _d_vector_i = clad::one_hot_vector(indepVarCount, {{0U|0UL|0ULL}});
+  // CHECK-NEXT:     clad::array<double> _d_vector_j = clad::one_hot_vector(indepVarCount, {{1U|1UL|1ULL}});
+  // CHECK-NEXT:     * _d_vector_output = clad::identity_matrix(* _d_vector_output.rows(), indepVarCount, {{2U|2UL|2ULL}});
+  // CHECK-NEXT:     volatile double &_t0 = this->x;
+  // CHECK-NEXT:     double _t1 = _t0 * i;
+  // CHECK-NEXT:     double _t2 = _t1 * i;
+  // CHECK-NEXT:     * _d_vector_output[0] = ((0 * i + _t0 * _d_vector_i) * i + _t1 * _d_vector_i) * j + _t2 * _d_vector_j;
+  // CHECK-NEXT:     output[0] = _t2 * j;
+  // CHECK-NEXT:     volatile double &_t3 = this->y;
+  // CHECK-NEXT:     double _t4 = _t3 * i;
+  // CHECK-NEXT:     double _t5 = _t4 * j;
+  // CHECK-NEXT:     * _d_vector_output[1] = ((0 * i + _t3 * _d_vector_i) * j + _t4 * _d_vector_j) * j + _t5 * _d_vector_j;
+  // CHECK-NEXT:     output[1] = _t5 * j;
   // CHECK-NEXT: }
 
 };
@@ -126,29 +106,21 @@ struct ExperimentConstVolatile {
     x = val;
   }
 
-  // CHECK: void operator_call_jac(double i, double j, double *output, double *jacobianMatrix) const volatile {
-  // CHECK-NEXT:     double _t0 = this->x * i;
-  // CHECK-NEXT:     double _t1 = output[0];
-  // CHECK-NEXT:     output[0] = this->x * i * i * j;
-  // CHECK-NEXT:     double _t2 = this->y * i;
-  // CHECK-NEXT:     double _t3 = output[1];
-  // CHECK-NEXT:     output[1] = this->y * i * j * j;
-  // CHECK-NEXT:     {
-  // CHECK-NEXT:         {
-  // CHECK-NEXT:             jacobianMatrix[{{2U|2UL|2ULL}}] += this->y * 1 * j * j;
-  // CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += _t2 * 1 * j;
-  // CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += _t2 * j * 1;
-  // CHECK-NEXT:         }
-  // CHECK-NEXT:         output[1] = _t3;
-  // CHECK-NEXT:     }
-  // CHECK-NEXT:     {
-  // CHECK-NEXT:         {
-  // CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += this->x * 1 * j * i;
-  // CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += _t0 * 1 * j;
-  // CHECK-NEXT:             jacobianMatrix[{{1U|1UL|1ULL}}] += _t0 * i * 1;
-  // CHECK-NEXT:         }
-  // CHECK-NEXT:         output[0] = _t1;
-  // CHECK-NEXT:     }
+  // CHECK: void operator_call_jac(double i, double j, double *output, clad::array_ref<clad::matrix<double> > _d_vector_output) const volatile {
+  // CHECK-NEXT:     unsigned long indepVarCount = * _d_vector_output.rows() + {{2U|2UL|2ULL}};
+  // CHECK-NEXT:     clad::array<double> _d_vector_i = clad::one_hot_vector(indepVarCount, {{0U|0UL|0ULL}});
+  // CHECK-NEXT:     clad::array<double> _d_vector_j = clad::one_hot_vector(indepVarCount, {{1U|1UL|1ULL}});
+  // CHECK-NEXT:     * _d_vector_output = clad::identity_matrix(* _d_vector_output.rows(), indepVarCount, {{2U|2UL|2ULL}});
+  // CHECK-NEXT:     volatile double &_t0 = this->x;
+  // CHECK-NEXT:     double _t1 = _t0 * i;
+  // CHECK-NEXT:     double _t2 = _t1 * i;
+  // CHECK-NEXT:     * _d_vector_output[0] = ((0 * i + _t0 * _d_vector_i) * i + _t1 * _d_vector_i) * j + _t2 * _d_vector_j;
+  // CHECK-NEXT:     output[0] = _t2 * j;
+  // CHECK-NEXT:     volatile double &_t3 = this->y;
+  // CHECK-NEXT:     double _t4 = _t3 * i;
+  // CHECK-NEXT:     double _t5 = _t4 * j;
+  // CHECK-NEXT:     * _d_vector_output[1] = ((0 * i + _t3 * _d_vector_i) * j + _t4 * _d_vector_j) * j + _t5 * _d_vector_j;
+  // CHECK-NEXT:     output[1] = _t5 * j;
   // CHECK-NEXT: }
 
 };
@@ -166,27 +138,21 @@ namespace outer {
         x = val;
       }
       
-  // CHECK: void operator_call_jac(double i, double j, double *output, double *jacobianMatrix) {
-  // CHECK-NEXT:     double _t0 = output[0];
-  // CHECK-NEXT:     output[0] = this->x * i * i * j;
-  // CHECK-NEXT:     double _t1 = output[1];
-  // CHECK-NEXT:     output[1] = this->y * i * j * j;
-  // CHECK-NEXT:     {
-  // CHECK-NEXT:         {
-  // CHECK-NEXT:             jacobianMatrix[{{2U|2UL|2ULL}}] += this->y * 1 * j * j;
-  // CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += this->y * i * 1 * j;
-  // CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += this->y * i * j * 1;
-  // CHECK-NEXT:         }
-  // CHECK-NEXT:         output[1] = _t1;
-  // CHECK-NEXT:     }
-  // CHECK-NEXT:     {
-  // CHECK-NEXT:         {
-  // CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += this->x * 1 * j * i;
-  // CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += this->x * i * 1 * j;
-  // CHECK-NEXT:             jacobianMatrix[{{1U|1UL|1ULL}}] += this->x * i * i * 1;
-  // CHECK-NEXT:         }
-  // CHECK-NEXT:         output[0] = _t0;
-  // CHECK-NEXT:     }
+  // CHECK: void operator_call_jac(double i, double j, double *output, clad::array_ref<clad::matrix<double> > _d_vector_output) {
+  // CHECK-NEXT:     unsigned long indepVarCount = * _d_vector_output.rows() + {{2U|2UL|2ULL}};
+  // CHECK-NEXT:     clad::array<double> _d_vector_i = clad::one_hot_vector(indepVarCount, {{0U|0UL|0ULL}});
+  // CHECK-NEXT:     clad::array<double> _d_vector_j = clad::one_hot_vector(indepVarCount, {{1U|1UL|1ULL}});
+  // CHECK-NEXT:     * _d_vector_output = clad::identity_matrix(* _d_vector_output.rows(), indepVarCount, {{2U|2UL|2ULL}});
+  // CHECK-NEXT:     double &_t0 = this->x;
+  // CHECK-NEXT:     double _t1 = _t0 * i;
+  // CHECK-NEXT:     double _t2 = _t1 * i;
+  // CHECK-NEXT:     * _d_vector_output[0] = ((0 * i + _t0 * _d_vector_i) * i + _t1 * _d_vector_i) * j + _t2 * _d_vector_j;
+  // CHECK-NEXT:     output[0] = _t2 * j;
+  // CHECK-NEXT:     double &_t3 = this->y;
+  // CHECK-NEXT:     double _t4 = _t3 * i;
+  // CHECK-NEXT:     double _t5 = _t4 * j;
+  // CHECK-NEXT:     * _d_vector_output[1] = ((0 * i + _t3 * _d_vector_i) * j + _t4 * _d_vector_j) * j + _t5 * _d_vector_j;
+  // CHECK-NEXT:     output[1] = _t5 * j;
   // CHECK-NEXT: }
 
     };
@@ -196,27 +162,17 @@ namespace outer {
       output[1] = i*j*j;
     };
 
-  // CHECK: inline void operator_call_jac(double i, double j, double *output, double *jacobianMatrix) const {
-  // CHECK-NEXT:     double _t0 = output[0];
-  // CHECK-NEXT:     output[0] = i * i * j;
-  // CHECK-NEXT:     double _t1 = output[1];
-  // CHECK-NEXT:     output[1] = i * j * j;
-  // CHECK-NEXT:     {
-  // CHECK-NEXT:         {
-  // CHECK-NEXT:             jacobianMatrix[{{2U|2UL|2ULL}}] += 1 * j * j;
-  // CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += i * 1 * j;
-  // CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += i * j * 1;
-  // CHECK-NEXT:         }
-  // CHECK-NEXT:         output[1] = _t1;
-  // CHECK-NEXT:     }
-  // CHECK-NEXT:     {
-  // CHECK-NEXT:         {
-  // CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += 1 * j * i;
-  // CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += i * 1 * j;
-  // CHECK-NEXT:             jacobianMatrix[{{1U|1UL|1ULL}}] += i * i * 1;
-  // CHECK-NEXT:         }
-  // CHECK-NEXT:         output[0] = _t0;
-  // CHECK-NEXT:     }
+  // CHECK: inline void operator_call_jac(double i, double j, double *output, clad::array_ref<clad::matrix<double> > _d_vector_output) const {
+  // CHECK-NEXT:     unsigned long indepVarCount = * _d_vector_output.rows() + {{2U|2UL|2ULL}};
+  // CHECK-NEXT:     clad::array<double> _d_vector_i = clad::one_hot_vector(indepVarCount, {{0U|0UL|0ULL}});
+  // CHECK-NEXT:     clad::array<double> _d_vector_j = clad::one_hot_vector(indepVarCount, {{1U|1UL|1ULL}});
+  // CHECK-NEXT:     * _d_vector_output = clad::identity_matrix(* _d_vector_output.rows(), indepVarCount, {{2U|2UL|2ULL}});
+  // CHECK-NEXT:     double _t0 = i * i;
+  // CHECK-NEXT:     * _d_vector_output[0] = (_d_vector_i * i + i * _d_vector_i) * j + _t0 * _d_vector_j;
+  // CHECK-NEXT:     output[0] = _t0 * j;
+  // CHECK-NEXT:     double _t1 = i * j;
+  // CHECK-NEXT:     * _d_vector_output[1] = (_d_vector_i * j + i * _d_vector_j) * j + _t1 * _d_vector_j;
+  // CHECK-NEXT:     output[1] = _t1 * j;
   // CHECK-NEXT: }
 
   }
@@ -227,21 +183,20 @@ namespace outer {
   auto d_##E##Ref = clad::jacobian(E);
 
 #define TEST(E)                                                                \
-  result[0] = result[1] = result[2] = result[3] = 0;                           \
   output[0] = output[1] = 0;                                                   \
-  d_##E.execute(7, 9, output, result);                                         \
-  printf("{%.2f, %.2f, %.2f, %.2f}, ", result[0], result[1], result[2],        \
-         result[3]);                                                           \
-  result[0] = result[1] = result[2] = result[3] = 0;                           \
+  d_##E.execute(7, 9, output, &result);                                        \
+  printf("{%.2f, %.2f, %.2f, %.2f}, ", result[0][0], result[0][1],             \
+                                       result[1][0], result[1][1]);            \
   output[0] = output[1] = 0;                                                   \
-  d_##E##Ref.execute(7, 9, output, result);                                    \
-  printf("{%.2f, %.2f, %.2f, %.2f}\n", result[0], result[1], result[2],        \
-         result[3]);
+  d_##E##Ref.execute(7, 9, output, &result);                                   \
+  printf("{%.2f, %.2f, %.2f, %.2f}, ", result[0][0], result[0][1],             \
+                                       result[1][0], result[1][1]);
 
 double x = 3;
 double y = 5;
 int main() {
-  double output[2], result[4];
+  double output[2];
+  clad::matrix<double> result(2, 2);
   Experiment E(3, 5);
   auto E_Again = E;
   const ExperimentConst E_Const(3, 5);
@@ -254,27 +209,17 @@ int main() {
     output[1] = i*j*j;
   };
 
-  // CHECK: inline void operator_call_jac(double i, double j, double *output, double *jacobianMatrix) const {
-  // CHECK-NEXT:     double _t0 = output[0];
-  // CHECK-NEXT:     output[0] = i * i * j;
-  // CHECK-NEXT:     double _t1 = output[1];
-  // CHECK-NEXT:     output[1] = i * j * j;
-  // CHECK-NEXT:     {
-  // CHECK-NEXT:         {
-  // CHECK-NEXT:             jacobianMatrix[{{2U|2UL|2ULL}}] += 1 * j * j;
-  // CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += i * 1 * j;
-  // CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += i * j * 1;
-  // CHECK-NEXT:         }
-  // CHECK-NEXT:         output[1] = _t1;
-  // CHECK-NEXT:     }
-  // CHECK-NEXT:     {
-  // CHECK-NEXT:         {
-  // CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += 1 * j * i;
-  // CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += i * 1 * j;
-  // CHECK-NEXT:             jacobianMatrix[{{1U|1UL|1ULL}}] += i * i * 1;
-  // CHECK-NEXT:         }
-  // CHECK-NEXT:         output[0] = _t0;
-  // CHECK-NEXT:     }
+  // CHECK-NEXT: inline void operator_call_jac(double i, double j, double *output, clad::array_ref<clad::matrix<double> > _d_vector_output) const {
+  // CHECK-NEXT:     unsigned long indepVarCount = * _d_vector_output.rows() + {{2U|2UL|2ULL}};
+  // CHECK-NEXT:     clad::array<double> _d_vector_i = clad::one_hot_vector(indepVarCount, {{0U|0UL|0ULL}});
+  // CHECK-NEXT:     clad::array<double> _d_vector_j = clad::one_hot_vector(indepVarCount, {{1U|1UL|1ULL}});
+  // CHECK-NEXT:     * _d_vector_output = clad::identity_matrix(* _d_vector_output.rows(), indepVarCount, {{2U|2UL|2ULL}});
+  // CHECK-NEXT:     double _t0 = i * i;
+  // CHECK-NEXT:     * _d_vector_output[0] = (_d_vector_i * i + i * _d_vector_i) * j + _t0 * _d_vector_j;
+  // CHECK-NEXT:     output[0] = _t0 * j;
+  // CHECK-NEXT:     double _t1 = i * j;
+  // CHECK-NEXT:     * _d_vector_output[1] = (_d_vector_i * j + i * _d_vector_j) * j + _t1 * _d_vector_j;
+  // CHECK-NEXT:     output[1] = _t1 * j;
   // CHECK-NEXT: }
 
   auto lambdaWithCapture = [&](double i, double jj, double *output) {
@@ -282,27 +227,19 @@ int main() {
     output[1] = y*i*jj*jj;
   };
 
-// CHECK: inline void operator_call_jac(double i, double jj, double *output, double *jacobianMatrix) const {
-// CHECK-NEXT:     double _t0 = output[0];
-// CHECK-NEXT:     output[0] = x * i * i * jj;
-// CHECK-NEXT:     double _t1 = output[1];
-// CHECK-NEXT:     output[1] = y * i * jj * jj;
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{2U|2UL|2ULL}}] += y * 1 * jj * jj;
-// CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += y * i * 1 * jj;
-// CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += y * i * jj * 1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         output[1] = _t1;
-// CHECK-NEXT:     }
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += x * 1 * jj * i;
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += x * i * 1 * jj;
-// CHECK-NEXT:             jacobianMatrix[{{1U|1UL|1ULL}}] += x * i * i * 1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         output[0] = _t0;
-// CHECK-NEXT:     }
+// CHECK: inline void operator_call_jac(double i, double jj, double *output, clad::array_ref<clad::matrix<double> > _d_vector_output) const {
+// CHECK-NEXT:     unsigned long indepVarCount = * _d_vector_output.rows() + {{2U|2UL|2ULL}};
+// CHECK-NEXT:     clad::array<double> _d_vector_i = clad::one_hot_vector(indepVarCount, {{0U|0UL|0ULL}});
+// CHECK-NEXT:     clad::array<double> _d_vector_jj = clad::one_hot_vector(indepVarCount, {{1U|1UL|1ULL}});
+// CHECK-NEXT:     * _d_vector_output = clad::identity_matrix(* _d_vector_output.rows(), indepVarCount, {{2U|2UL|2ULL}});
+// CHECK-NEXT:     double _t0 = x * i;
+// CHECK-NEXT:     double _t1 = _t0 * i;
+// CHECK-NEXT:     * _d_vector_output[0] = ((0. * i + x * _d_vector_i) * i + _t0 * _d_vector_i) * jj + _t1 * _d_vector_jj;
+// CHECK-NEXT:     output[0] = _t1 * jj;
+// CHECK-NEXT:     double _t2 = y * i;
+// CHECK-NEXT:     double _t3 = _t2 * jj;
+// CHECK-NEXT:     * _d_vector_output[1] = ((0. * i + y * _d_vector_i) * jj + _t2 * _d_vector_jj) * jj + _t3 * _d_vector_jj;
+// CHECK-NEXT:     output[1] = _t3 * jj;
 // CHECK-NEXT: }
 
   auto lambdaNNS = outer::inner::lambdaNNS;
diff --git a/test/Jacobian/Jacobian.C b/test/Jacobian/Jacobian.C
index 4d9537022..9d8eeb298 100644
--- a/test/Jacobian/Jacobian.C
+++ b/test/Jacobian/Jacobian.C
@@ -12,43 +12,22 @@ void f_1(double a, double b, double c, double output[]) {
   output[2] = c * c * 10 - a * a;
 }
 
-void f_1_jac(double a, double b, double c, double output[], double *_result);
-
-// CHECK: void f_1_jac(double a, double b, double c, double output[], double *jacobianMatrix) {
-// CHECK-NEXT:     double _t0 = output[0];
-// CHECK-NEXT:     output[0] = a * a * a;
-// CHECK-NEXT:     double _t1 = output[1];
-// CHECK-NEXT:     output[1] = a * a * a + b * b * b;
-// CHECK-NEXT:     double _t2 = output[2];
-// CHECK-NEXT:     output[2] = c * c * 10 - a * a;
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{8U|8UL|8ULL}}] += 1 * 10 * c;
-// CHECK-NEXT:             jacobianMatrix[{{8U|8UL|8ULL}}] += c * 1 * 10;
-// CHECK-NEXT:             jacobianMatrix[{{6U|6UL|6ULL}}] += -1 * a;
-// CHECK-NEXT:             jacobianMatrix[{{6U|6UL|6ULL}}] += a * -1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         output[2] = _t2;
-// CHECK-NEXT:     }
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += 1 * a * a;
-// CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += a * 1 * a;
-// CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += a * a * 1;
-// CHECK-NEXT:             jacobianMatrix[{{4U|4UL|4ULL}}] += 1 * b * b;
-// CHECK-NEXT:             jacobianMatrix[{{4U|4UL|4ULL}}] += b * 1 * b;
-// CHECK-NEXT:             jacobianMatrix[{{4U|4UL|4ULL}}] += b * b * 1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         output[1] = _t1;
-// CHECK-NEXT:     }
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += 1 * a * a;
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += a * 1 * a;
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += a * a * 1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         output[0] = _t0;
-// CHECK-NEXT:     }
+// CHECK: void f_1_jac(double a, double b, double c, double output[], clad::array_ref<clad::matrix<double> > _d_vector_output) {
+// CHECK-NEXT:     unsigned long indepVarCount = * _d_vector_output.rows() + {{3U|3UL|3ULL}};
+// CHECK-NEXT:     clad::array<double> _d_vector_a = clad::one_hot_vector(indepVarCount, {{0U|0UL|0ULL}});
+// CHECK-NEXT:     clad::array<double> _d_vector_b = clad::one_hot_vector(indepVarCount, {{1U|1UL|1ULL}});
+// CHECK-NEXT:     clad::array<double> _d_vector_c = clad::one_hot_vector(indepVarCount, {{2U|2UL|2ULL}});
+// CHECK-NEXT:     * _d_vector_output = clad::identity_matrix(* _d_vector_output.rows(), indepVarCount, {{3U|3UL|3ULL}});
+// CHECK-NEXT:     double _t0 = a * a;
+// CHECK-NEXT:     * _d_vector_output[0] = (_d_vector_a * a + a * _d_vector_a) * a + _t0 * _d_vector_a;
+// CHECK-NEXT:     output[0] = _t0 * a;
+// CHECK-NEXT:     double _t1 = a * a;
+// CHECK-NEXT:     double _t2 = b * b;
+// CHECK-NEXT:     * _d_vector_output[1] = (_d_vector_a * a + a * _d_vector_a) * a + _t1 * _d_vector_a + (_d_vector_b * b + b * _d_vector_b) * b + _t2 * _d_vector_b;
+// CHECK-NEXT:     output[1] = _t1 * a + _t2 * b;
+// CHECK-NEXT:     double _t3 = c * c;
+// CHECK-NEXT:     * _d_vector_output[2] = (_d_vector_c * c + c * _d_vector_c) * 10 + _t3 * 0 - (_d_vector_a * a + a * _d_vector_a);
+// CHECK-NEXT:     output[2] = _t3 * 10 - a * a;
 // CHECK-NEXT: }
 
 void f_3(double x, double y, double z, double *_result) {
@@ -59,47 +38,29 @@ void f_3(double x, double y, double z, double *_result) {
   _result[2] = sin(z) * constant;
 }
 
-void f_3_jac(double x, double y, double z, double *_result, double *jacobianMatrix);
-
-// CHECK: void f_3_jac(double x, double y, double z, double *_result, double *jacobianMatrix) {
-// CHECK-NEXT:     double _d_constant = 0.;
+// CHECK: void f_3_jac(double x, double y, double z, double *_result, clad::array_ref<clad::matrix<double> > _d_vector__result) {
+// CHECK-NEXT:     unsigned long indepVarCount = * _d_vector__result.rows() + {{3U|3UL|3ULL}};
+// CHECK-NEXT:     clad::array<double> _d_vector_x = clad::one_hot_vector(indepVarCount, {{0U|0UL|0ULL}});
+// CHECK-NEXT:     clad::array<double> _d_vector_y = clad::one_hot_vector(indepVarCount, {{1U|1UL|1ULL}});
+// CHECK-NEXT:     clad::array<double> _d_vector_z = clad::one_hot_vector(indepVarCount, {{2U|2UL|2ULL}});
+// CHECK-NEXT:     * _d_vector__result = clad::identity_matrix(* _d_vector__result.rows(), indepVarCount, {{3U|3UL|3ULL}});
+// CHECK-NEXT:     clad::array<double> _d_vector_constant(clad::array<double>(indepVarCount, 0));
 // CHECK-NEXT:     double constant = 42;
-// CHECK-NEXT:     double _t0 = sin(x);
-// CHECK-NEXT:     double _t1 = _result[0];
-// CHECK-NEXT:     _result[0] = sin(x) * constant;
-// CHECK-NEXT:     double _t2 = sin(y);
-// CHECK-NEXT:     double _t3 = _result[1];
-// CHECK-NEXT:     _result[1] = sin(y) * constant;
-// CHECK-NEXT:     double _t4 = sin(z);
-// CHECK-NEXT:     double _t5 = _result[2];
-// CHECK-NEXT:     _result[2] = sin(z) * constant;
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             double _r2 = 0.;
-// CHECK-NEXT:             _r2 += 1 * constant * clad::custom_derivatives::sin_pushforward(z, 1.).pushforward;
-// CHECK-NEXT:             jacobianMatrix[{{8U|8UL|8ULL}}] += _r2;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         _result[2] = _t5;
-// CHECK-NEXT:     }
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             double _r1 = 0.;
-// CHECK-NEXT:             _r1 += 1 * constant * clad::custom_derivatives::sin_pushforward(y, 1.).pushforward;
-// CHECK-NEXT:             jacobianMatrix[{{4U|4UL|4ULL}}] += _r1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         _result[1] = _t3;
-// CHECK-NEXT:     }
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             double _r0 = 0.;
-// CHECK-NEXT:             _r0 += 1 * constant * clad::custom_derivatives::sin_pushforward(x, 1.).pushforward;
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += _r0;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         _result[0] = _t1;
-// CHECK-NEXT:     }
+// CHECK-NEXT:     {{.*}} _t0 = clad::custom_derivatives::sin_pushforward(x, _d_vector_x);
+// CHECK-NEXT:     double &_t1 = _t0.value;
+// CHECK-NEXT:     * _d_vector__result[0] = _t0.pushforward * constant + _t1 * _d_vector_constant;
+// CHECK-NEXT:     _result[0] = _t1 * constant;
+// CHECK-NEXT:     {{.*}} _t2 = clad::custom_derivatives::sin_pushforward(y, _d_vector_y);
+// CHECK-NEXT:     double &_t3 = _t2.value;
+// CHECK-NEXT:     * _d_vector__result[1] = _t2.pushforward * constant + _t3 * _d_vector_constant;
+// CHECK-NEXT:     _result[1] = _t3 * constant;
+// CHECK-NEXT:     {{.*}} _t4 = clad::custom_derivatives::sin_pushforward(z, _d_vector_z);
+// CHECK-NEXT:     double &_t5 = _t4.value;
+// CHECK-NEXT:     * _d_vector__result[2] = _t4.pushforward * constant + _t5 * _d_vector_constant;
+// CHECK-NEXT:     _result[2] = _t5 * constant;
 // CHECK-NEXT: }
 double multiply(double x, double y) { return x * y; }
-//CHECK: void multiply_pullback(double x, double y, double _d_y0, double *_d_x, double *_d_y);
+// CHECK: clad::ValueAndPushforward<double, clad::array<double> > multiply_vector_pushforward(double x, double y, clad::array<double> _d_x, clad::array<double> _d_y);
 
 void f_4(double x, double y, double z, double *_result) {
   double constant = 42;
@@ -109,84 +70,43 @@ void f_4(double x, double y, double z, double *_result) {
   _result[2] = multiply(z, x) * constant;
 }
 
-void f_4_jac(double x, double y, double z, double *_result, double *jacobianMatrix);
-// CHECK: void f_4_jac(double x, double y, double z, double *_result, double *jacobianMatrix) {
-// CHECK-NEXT:     double _d_constant = 0.;
+// CHECK: void f_4_jac(double x, double y, double z, double *_result, clad::array_ref<clad::matrix<double> > _d_vector__result) {
+// CHECK-NEXT:     unsigned long indepVarCount = * _d_vector__result.rows() + {{3U|3UL|3ULL}};
+// CHECK-NEXT:     clad::array<double> _d_vector_x = clad::one_hot_vector(indepVarCount, {{0U|0UL|0ULL}});
+// CHECK-NEXT:     clad::array<double> _d_vector_y = clad::one_hot_vector(indepVarCount, {{1U|1UL|1ULL}});
+// CHECK-NEXT:     clad::array<double> _d_vector_z = clad::one_hot_vector(indepVarCount, {{2U|2UL|2ULL}});
+// CHECK-NEXT:     * _d_vector__result = clad::identity_matrix(* _d_vector__result.rows(), indepVarCount, {{3U|3UL|3ULL}});
+// CHECK-NEXT:     clad::array<double> _d_vector_constant(clad::array<double>(indepVarCount, 0));
 // CHECK-NEXT:     double constant = 42;
-// CHECK-NEXT:     double _t0 = multiply(x, y);
-// CHECK-NEXT:     double _t1 = _result[0];
-// CHECK-NEXT:     _result[0] = multiply(x, y) * constant;
-// CHECK-NEXT:     double _t2 = multiply(y, z);
-// CHECK-NEXT:     double _t3 = _result[1];
-// CHECK-NEXT:     _result[1] = multiply(y, z) * constant;
-// CHECK-NEXT:     double _t4 = multiply(z, x);
-// CHECK-NEXT:     double _t5 = _result[2];
-// CHECK-NEXT:     _result[2] = multiply(z, x) * constant;
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             double _r4 = 0.;
-// CHECK-NEXT:             double _r5 = 0.;
-// CHECK-NEXT:             multiply_pullback(z, x, 1 * constant, &_r4, &_r5);
-// CHECK-NEXT:             jacobianMatrix[{{8U|8UL|8ULL}}] += _r4;
-// CHECK-NEXT:             jacobianMatrix[{{6U|6UL|6ULL}}] += _r5;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         _result[2] = _t5;
-// CHECK-NEXT:     }
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             double _r2 = 0.;
-// CHECK-NEXT:             double _r3 = 0.;
-// CHECK-NEXT:             multiply_pullback(y, z, 1 * constant, &_r2, &_r3);
-// CHECK-NEXT:             jacobianMatrix[{{4U|4UL|4ULL}}] += _r2;
-// CHECK-NEXT:             jacobianMatrix[{{5U|5UL|5ULL}}] += _r3;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         _result[1] = _t3;
-// CHECK-NEXT:     }
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             double _r0 = 0.;
-// CHECK-NEXT:             double _r1 = 0.;
-// CHECK-NEXT:             multiply_pullback(x, y, 1 * constant, &_r0, &_r1);
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += _r0;
-// CHECK-NEXT:             jacobianMatrix[{{1U|1UL|1ULL}}] += _r1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         _result[0] = _t1;
-// CHECK-NEXT:     }
+// CHECK-NEXT:     clad::ValueAndPushforward<double, clad::array<double> > _t0 = multiply_vector_pushforward(x, y, _d_vector_x, _d_vector_y);
+// CHECK-NEXT:     double &_t1 = _t0.value;
+// CHECK-NEXT:     * _d_vector__result[0] = _t0.pushforward * constant + _t1 * _d_vector_constant;
+// CHECK-NEXT:     _result[0] = _t1 * constant;
+// CHECK-NEXT:     clad::ValueAndPushforward<double, clad::array<double> > _t2 = multiply_vector_pushforward(y, z, _d_vector_y, _d_vector_z);
+// CHECK-NEXT:     double &_t3 = _t2.value;
+// CHECK-NEXT:     * _d_vector__result[1] = _t2.pushforward * constant + _t3 * _d_vector_constant;
+// CHECK-NEXT:     _result[1] = _t3 * constant;
+// CHECK-NEXT:     clad::ValueAndPushforward<double, clad::array<double> > _t4 = multiply_vector_pushforward(z, x, _d_vector_z, _d_vector_x);
+// CHECK-NEXT:     double &_t5 = _t4.value;
+// CHECK-NEXT:     * _d_vector__result[2] = _t4.pushforward * constant + _t5 * _d_vector_constant;
+// CHECK-NEXT:     _result[2] = _t5 * constant;
 // CHECK-NEXT: }
 
-void f_1_jac_0(double a, double b, double c, double output[], double *jacobianMatrix);
-// CHECK: void f_1_jac_0(double a, double b, double c, double output[], double *jacobianMatrix) {
-// CHECK-NEXT:     double _d_b = 0.;
-// CHECK-NEXT:     double _d_c = 0.;
-// CHECK-NEXT:     double _t0 = output[0];
-// CHECK-NEXT:     output[0] = a * a * a;
-// CHECK-NEXT:     double _t1 = output[1];
-// CHECK-NEXT:     output[1] = a * a * a + b * b * b;
-// CHECK-NEXT:     double _t2 = output[2];
-// CHECK-NEXT:     output[2] = c * c * 10 - a * a;
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{2U|2UL|2ULL}}] += -1 * a;
-// CHECK-NEXT:             jacobianMatrix[{{2U|2UL|2ULL}}] += a * -1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         output[2] = _t2;
-// CHECK-NEXT:     }
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{1U|1UL|1ULL}}] += 1 * a * a;
-// CHECK-NEXT:             jacobianMatrix[{{1U|1UL|1ULL}}] += a * 1 * a;
-// CHECK-NEXT:             jacobianMatrix[{{1U|1UL|1ULL}}] += a * a * 1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         output[1] = _t1;
-// CHECK-NEXT:     }
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += 1 * a * a;
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += a * 1 * a;
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += a * a * 1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         output[0] = _t0;
-// CHECK-NEXT:     }
+// CHECK: void f_1_jac_0(double a, double b, double c, double output[], clad::array_ref<clad::matrix<double> > _d_vector_output) {
+// CHECK-NEXT:     unsigned long indepVarCount = * _d_vector_output.rows() + {{3U|3UL|3ULL}};
+// CHECK-NEXT:     clad::array<double> _d_vector_a = clad::one_hot_vector(indepVarCount, {{0U|0UL|0ULL}});
+// CHECK-NEXT:     clad::array<double> _d_vector_b = clad::zero_vector(indepVarCount);
+// CHECK-NEXT:     clad::array<double> _d_vector_c = clad::zero_vector(indepVarCount);
+// CHECK-NEXT:     double _t0 = a * a;
+// CHECK-NEXT:     * _d_vector_output[0] = (_d_vector_a * a + a * _d_vector_a) * a + _t0 * _d_vector_a;
+// CHECK-NEXT:     output[0] = _t0 * a;
+// CHECK-NEXT:     double _t1 = a * a;
+// CHECK-NEXT:     double _t2 = b * b;
+// CHECK-NEXT:     * _d_vector_output[1] = (_d_vector_a * a + a * _d_vector_a) * a + _t1 * _d_vector_a + (_d_vector_b * b + b * _d_vector_b) * b + _t2 * _d_vector_b;
+// CHECK-NEXT:     output[1] = _t1 * a + _t2 * b;
+// CHECK-NEXT:     double _t3 = c * c;
+// CHECK-NEXT:     * _d_vector_output[2] = (_d_vector_c * c + c * _d_vector_c) * 10 + _t3 * 0 - (_d_vector_a * a + a * _d_vector_a);
+// CHECK-NEXT:     output[2] = _t3 * 10 - a * a;
 // CHECK-NEXT: }
 
 void f_5(float a, double output[]){
@@ -194,51 +114,38 @@ void f_5(float a, double output[]){
   output[0]=a*a;  
 }
 
-// CHECK: void f_5_jac(float a, double output[], double *jacobianMatrix) {
-// CHECK-NEXT:     double _t0 = output[1];
+// CHECK: void f_5_jac(float a, double output[], clad::array_ref<clad::matrix<double> > _d_vector_output) {
+// CHECK-NEXT:     unsigned long indepVarCount = * _d_vector_output.rows() + {{1U|1UL|1ULL}};
+// CHECK-NEXT:     clad::array<float> _d_vector_a = clad::one_hot_vector(indepVarCount, {{0U|0UL|0ULL}});
+// CHECK-NEXT:     * _d_vector_output = clad::identity_matrix(* _d_vector_output.rows(), indepVarCount, {{1U|1UL|1ULL}});
+// CHECK-NEXT:     * _d_vector_output[1] = _d_vector_a;
 // CHECK-NEXT:     output[1] = a;
-// CHECK-NEXT:     double _t1 = output[0];
+// CHECK-NEXT:     * _d_vector_output[0] = _d_vector_a * a + a * _d_vector_a;
 // CHECK-NEXT:     output[0] = a * a;
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += 1 * a;
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += a * 1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         output[0] = _t1;
-// CHECK-NEXT:     }
-// CHECK-NEXT:     {
-// CHECK-NEXT:         jacobianMatrix[{{1U|1UL|1ULL}}] += 1;
-// CHECK-NEXT:         output[1] = _t0;
-// CHECK-NEXT:     }
 // CHECK-NEXT: }
 
 
 #define TEST(F, x, y, z) { \
-  result[0] = 0; result[1] = 0; result[2] = 0;\
-  result[3] = 0; result[4] = 0; result[5] = 0;\
-  result[6] = 0; result[7] = 0; result[8] = 0;\
   outputarr[0] = 0; outputarr[1] = 1; outputarr[2] = 0;\
   auto j = clad::jacobian(F);\
-  j.execute(x, y, z, outputarr, result);\
+  j.execute(x, y, z, outputarr, &result);\
   printf("Result is = {%.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f}\n",\
-  result[0], result[1], result[2],\
-  result[3], result[4], result[5],\
-  result[6], result[7], result[8]);\
-  F##_jac(x, y, z, outputarr, result);\
+  result[0][0], result[0][1], result[0][2],\
+  result[1][0], result[1][1], result[1][2],\
+  result[2][0], result[2][1], result[2][2]);\
 }
 
 #define TEST_F_1_SINGLE_PARAM(x, y, z) { \
-  result[0] = 0; result[1] = 0; result[2] = 0;\
   outputarr[0] = 0; outputarr[1] = 1; outputarr[2] = 0;\
   auto j = clad::jacobian(f_1,"a");\
-  j.execute(x, y, z, outputarr, result);\
+  j.execute(x, y, z, outputarr, &result);\
   printf("Result is = {%.2f, %.2f, %.2f}\n",\
-  result[0], result[1], result[2]);\
+  result[0][0], result[1][0], result[2][0]);\
 }
 
 
 int main() {
-  double result[10];
+  clad::matrix<double> result (3, 3);
   double outputarr[9];
   TEST(f_1, 1, 2, 3); // CHECK-EXEC: Result is = {3.00, 0.00, 0.00, 3.00, 12.00, 0.00, -2.00, 0.00, 60.00}
   TEST(f_3, 1, 2, 3); // CHECK-EXEC: Result is = {22.69, 0.00, 0.00, 0.00, -17.48, 0.00, 0.00, 0.00, -41.58}
@@ -247,13 +154,11 @@ int main() {
 
   auto df5 = clad::jacobian(f_5);
   result[0] = 0; result[1] = 0;
-  df5.execute(3, outputarr, result);
-  printf("Result is = {%.2f, %.2f}", result[0], result[1]); // CHECK-EXEC: Result is = {6.00, 1.00}
+  df5.execute(3, outputarr, &result);
+  printf("Result is = {%.2f, %.2f}", result[0][0], result[1][0]); // CHECK-EXEC: Result is = {6.00, 1.00}
 }
 
-//CHECK: void multiply_pullback(double x, double y, double _d_y0, double *_d_x, double *_d_y) {
-//CHECK-NEXT:    {
-//CHECK-NEXT:        *_d_x += _d_y0 * y;
-//CHECK-NEXT:        *_d_y += x * _d_y0;
-//CHECK-NEXT:    }
-//CHECK-NEXT:}
+// CHECK: clad::ValueAndPushforward<double, clad::array<double> > multiply_vector_pushforward(double x, double y, clad::array<double> _d_x, clad::array<double> _d_y) {
+// CHECK-NEXT:     unsigned long indepVarCount = _d_y.size();
+// CHECK-NEXT:     return {x * y, _d_x * y + x * _d_y};
+// CHECK-NEXT: }
diff --git a/test/Jacobian/Pointers.C b/test/Jacobian/Pointers.C
index 01bdcecc4..ec2154ee7 100644
--- a/test/Jacobian/Pointers.C
+++ b/test/Jacobian/Pointers.C
@@ -10,31 +10,27 @@ void nonMemFn(double i, double j, double* out) {
   out[1] = j;
 }
 
-// CHECK: void nonMemFn_jac(double i, double j, double *out, double *jacobianMatrix) {
-// CHECK-NEXT:     double _t0 = out[0];
+// CHECK: void nonMemFn_jac(double i, double j, double *out, clad::array_ref<clad::matrix<double> > _d_vector_out) {
+// CHECK-NEXT:     unsigned long indepVarCount = * _d_vector_out.rows() + {{2U|2UL|2ULL}};
+// CHECK-NEXT:     clad::array<double> _d_vector_i = clad::one_hot_vector(indepVarCount, {{0U|0UL|0ULL}});
+// CHECK-NEXT:     clad::array<double> _d_vector_j = clad::one_hot_vector(indepVarCount, {{1U|1UL|1ULL}});
+// CHECK-NEXT:     * _d_vector_out = clad::identity_matrix(* _d_vector_out.rows(), indepVarCount, {{2U|2UL|2ULL}});
+// CHECK-NEXT:     * _d_vector_out[0] = _d_vector_i;
 // CHECK-NEXT:     out[0] = i;
-// CHECK-NEXT:     double _t1 = out[1];
+// CHECK-NEXT:     * _d_vector_out[1] = _d_vector_j;
 // CHECK-NEXT:     out[1] = j;
-// CHECK-NEXT:     {
-// CHECK-NEXT:         jacobianMatrix[{{3U|3UL|3ULL}}] += 1;
-// CHECK-NEXT:         out[1] = _t1;
-// CHECK-NEXT:     }
-// CHECK-NEXT:     {
-// CHECK-NEXT:         jacobianMatrix[{{0U|0UL|0ULL}}] += 1;
-// CHECK-NEXT:         out[0] = _t0;
-// CHECK-NEXT:     }
 // CHECK-NEXT: }
 
 #define NON_MEM_FN_TEST(var)\
-res[0]=res[1]=res[2]=res[3]=0;\
-var.execute(5, 7, out, res);\
-printf("{%.2f %.2f %.2f %.2f}\n", res[0], res[1], res[2], res[3]);
+var.execute(5, 7, out, &res);\
+printf("{%.2f %.2f %.2f %.2f}\n", res[0][0], res[0][1],\
+                                  res[1][0], res[1][1]);
 
 int main() {
   auto nonMemFnPtr = &nonMemFn;
   auto nonMemFnPtrToPtr = &nonMemFnPtr;
 
-  double res[4];
+  clad::matrix<double> res(2, 2);
   double out[2];
   auto d_nonMemFn = clad::jacobian(nonMemFn);
   auto d_nonMemFnPar = clad::jacobian((nonMemFn));
diff --git a/test/Jacobian/TemplateFunctors.C b/test/Jacobian/TemplateFunctors.C
index 9762d355e..f5141e0ad 100644
--- a/test/Jacobian/TemplateFunctors.C
+++ b/test/Jacobian/TemplateFunctors.C
@@ -15,25 +15,24 @@ template <typename T> struct Experiment {
   void setX(T val) { x = val; }
 };
 
-// CHECK: void operator_call_jac(double i, double j, double *output, double *jacobianMatrix) {
-// CHECK-NEXT:     double _t0 = output[0];
-// CHECK-NEXT:     output[0] = this->x * this->y * i * j;
-// CHECK-NEXT:     double _t1 = output[1];
-// CHECK-NEXT:     output[1] = 2 * this->x * this->y * i * j;
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{2U|2UL|2ULL}}] += 2 * this->x * this->y * 1 * j;
-// CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += 2 * this->x * this->y * i * 1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         output[1] = _t1;
-// CHECK-NEXT:     }
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += this->x * this->y * 1 * j;
-// CHECK-NEXT:             jacobianMatrix[{{1U|1UL|1ULL}}] += this->x * this->y * i * 1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         output[0] = _t0;
-// CHECK-NEXT:     }
+// CHECK: void operator_call_jac(double i, double j, double *output, clad::array_ref<clad::matrix<double> > _d_vector_output) {
+// CHECK-NEXT:     unsigned long indepVarCount = * _d_vector_output.rows() + {{2U|2UL|2ULL}};
+// CHECK-NEXT:     clad::array<double> _d_vector_i = clad::one_hot_vector(indepVarCount, {{0U|0UL|0ULL}});
+// CHECK-NEXT:     clad::array<double> _d_vector_j = clad::one_hot_vector(indepVarCount, {{1U|1UL|1ULL}});
+// CHECK-NEXT:     * _d_vector_output = clad::identity_matrix(* _d_vector_output.rows(), indepVarCount, {{2U|2UL|2ULL}});
+// CHECK-NEXT:     double &_t0 = this->x;
+// CHECK-NEXT:     double &_t1 = this->y;
+// CHECK-NEXT:     double _t2 = _t0 * _t1;
+// CHECK-NEXT:     double _t3 = _t2 * i;
+// CHECK-NEXT:     * _d_vector_output[0] = ((0 * _t1 + _t0 * 0) * i + _t2 * _d_vector_i) * j + _t3 * _d_vector_j;
+// CHECK-NEXT:     output[0] = _t3 * j;
+// CHECK-NEXT:     double &_t4 = this->x;
+// CHECK-NEXT:     double _t5 = 2 * _t4;
+// CHECK-NEXT:     double &_t6 = this->y;
+// CHECK-NEXT:     double _t7 = _t5 * _t6;
+// CHECK-NEXT:     double _t8 = _t7 * i;
+// CHECK-NEXT:     * _d_vector_output[1] = (((0 * _t4 + 2 * 0) * _t6 + _t5 * 0) * i + _t7 * _d_vector_i) * j + _t8 * _d_vector_j;
+// CHECK-NEXT:     output[1] = _t8 * j;
 // CHECK-NEXT: }
 
 template <> struct Experiment<long double> {
@@ -46,27 +45,26 @@ template <> struct Experiment<long double> {
   void setX(long double val) { x = val; }
 };
 
-// CHECK: void operator_call_jac(long double i, long double j, long double *output, long double *jacobianMatrix) {
-// CHECK-NEXT:     long double _t0 = output[0];
-// CHECK-NEXT:     output[0] = this->x * this->y * i * i * j;
-// CHECK-NEXT:     long double _t1 = output[1];
-// CHECK-NEXT:     output[1] = 2 * this->x * this->y * i * i * j;
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{2U|2UL|2ULL}}] += 2 * this->x * this->y * 1 * j * i;
-// CHECK-NEXT:             jacobianMatrix[{{2U|2UL|2ULL}}] += 2 * this->x * this->y * i * 1 * j;
-// CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += 2 * this->x * this->y * i * i * 1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         output[1] = _t1;
-// CHECK-NEXT:     }
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += this->x * this->y * 1 * j * i;
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += this->x * this->y * i * 1 * j;
-// CHECK-NEXT:             jacobianMatrix[{{1U|1UL|1ULL}}] += this->x * this->y * i * i * 1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         output[0] = _t0;
-// CHECK-NEXT:     }
+// CHECK: void operator_call_jac(long double i, long double j, long double *output, clad::array_ref<clad::matrix<long double> > _d_vector_output) {
+// CHECK-NEXT:     unsigned long indepVarCount = * _d_vector_output.rows() + {{2U|2UL|2ULL}};
+// CHECK-NEXT:     clad::array<long double> _d_vector_i = clad::one_hot_vector(indepVarCount, {{0U|0UL|0ULL}});
+// CHECK-NEXT:     clad::array<long double> _d_vector_j = clad::one_hot_vector(indepVarCount, {{1U|1UL|1ULL}});
+// CHECK-NEXT:     * _d_vector_output = clad::identity_matrix(* _d_vector_output.rows(), indepVarCount, {{2U|2UL|2ULL}});
+// CHECK-NEXT:     long double &_t0 = this->x;
+// CHECK-NEXT:     long double &_t1 = this->y;
+// CHECK-NEXT:     long double _t2 = _t0 * _t1;
+// CHECK-NEXT:     long double _t3 = _t2 * i;
+// CHECK-NEXT:     long double _t4 = _t3 * i;
+// CHECK-NEXT:     * _d_vector_output[0] = (((0 * _t1 + _t0 * 0) * i + _t2 * _d_vector_i) * i + _t3 * _d_vector_i) * j + _t4 * _d_vector_j;
+// CHECK-NEXT:     output[0] = _t4 * j;
+// CHECK-NEXT:     long double &_t5 = this->x;
+// CHECK-NEXT:     long double _t6 = 2 * _t5;
+// CHECK-NEXT:     long double &_t7 = this->y;
+// CHECK-NEXT:     long double _t8 = _t6 * _t7;
+// CHECK-NEXT:     long double _t9 = _t8 * i;
+// CHECK-NEXT:     long double _t10 = _t9 * i;
+// CHECK-NEXT:     * _d_vector_output[1] = ((((0 * _t5 + 2 * 0) * _t7 + _t6 * 0) * i + _t8 * _d_vector_i) * i + _t9 * _d_vector_i) * j + _t10 * _d_vector_j;
+// CHECK-NEXT:     output[1] = _t10 * j;
 // CHECK-NEXT: }
 
 #define INIT(E)                                                                \
@@ -74,32 +72,30 @@ template <> struct Experiment<long double> {
   auto d_##E##Ref = clad::jacobian(E);
 
 #define TEST_DOUBLE(E, ...)                                                    \
-  result[0] = result[1] = result[2] = result[3] = 0;                           \
   output[0] = output[1] = 0;                                                   \
-  d_##E.execute(__VA_ARGS__, output, result);                                  \
-  printf("{%.2f, %.2f, %.2f, %.2f} ", result[0], result[1], result[2],        \
-         result[3]);                                                           \
-  result[0] = result[1] = result[2] = result[3] = 0;                           \
+  d_##E.execute(__VA_ARGS__, output, &result);                                 \
+  printf("{%.2f, %.2f, %.2f, %.2f} ", result[0][0], result[0][1],              \
+                                      result[1][0], result[1][1]);             \
   output[0] = output[1] = 0;                                                   \
-  d_##E##Ref.execute(__VA_ARGS__, output, result);                             \
-  printf("{%.2f, %.2f, %.2f, %.2f}\n", result[0], result[1], result[2],        \
-         result[3]);
+  d_##E##Ref.execute(__VA_ARGS__, output, &result);                            \
+  printf("{%.2f, %.2f, %.2f, %.2f} ", result[0][0], result[0][1],              \
+                                      result[1][0], result[1][1]);             
 
 #define TEST_LONG_DOUBLE(E, ...)                                               \
-  result_ld[0] = result_ld[1] = result_ld[2] = result_ld[3] = 0;               \
   output_ld[0] = output_ld[1] = 0;                                             \
-  d_##E.execute(__VA_ARGS__, output_ld, result_ld);                            \
-  printf("{%.2Lf, %.2Lf, %.2Lf, %.2Lf} ", result_ld[0], result_ld[1],         \
-         result_ld[2], result_ld[3]);                                          \
-  result_ld[0] = result_ld[1] = result_ld[2] = result_ld[3] = 0;               \
+  d_##E.execute(__VA_ARGS__, output_ld, &result_ld);                           \
+  printf("{%.2Lf, %.2Lf, %.2Lf, %.2Lf} ", result_ld[0][0], result_ld[0][1],    \
+                                      result_ld[1][0], result_ld[1][1]);       \
   output_ld[0] = output_ld[1] = 0;                                             \
-  d_##E##Ref.execute(__VA_ARGS__, output_ld, result_ld);                       \
-  printf("{%.2Lf, %.2Lf, %.2Lf, %.2Lf}\n", result_ld[0], result_ld[1],         \
-         result_ld[2], result_ld[3]);
+  d_##E##Ref.execute(__VA_ARGS__, output_ld, &result_ld);                      \
+  printf("{%.2Lf, %.2Lf, %.2Lf, %.2Lf} ", result_ld[0][0], result_ld[0][1],    \
+                                      result_ld[1][0], result_ld[1][1]);             
 
 int main() {
-  double result[4], output[2];
-  long double result_ld[4], output_ld[2];
+  double output[2];
+  clad::matrix<double> result(2, 2);
+  long double output_ld[2];
+  clad::matrix<long double> result_ld(2, 2);
   Experiment<double> E(3, 5);
   Experiment<long double> E_ld(3, 5);
 
diff --git a/test/Jacobian/constexprTest.C b/test/Jacobian/constexprTest.C
index 1d34cd196..25a7210da 100644
--- a/test/Jacobian/constexprTest.C
+++ b/test/Jacobian/constexprTest.C
@@ -8,9 +8,8 @@
 #include "../TestUtils.h"
 
   double result[3] = {0};
-  double jacobianou[6] = {0};
   double result1[3] = {0};
-  double jacobianou1[9] = {0};
+  clad::matrix<double> jacobian(3, 2);
 
 constexpr void fn_mul(double i, double j, double *res) {
    res[0] = i*i;
@@ -18,34 +17,17 @@ constexpr void fn_mul(double i, double j, double *res) {
    res[2] = i*j;
 }
 
-// CHECK: void fn_mul_jac(double i, double j, double *res, double *jacobianMatrix) {
-// CHECK-NEXT:     double _t0 = res[0];
+// CHECK: constexpr void fn_mul_jac(double i, double j, double *res, clad::array_ref<clad::matrix<double> > _d_vector_res) {
+// CHECK-NEXT:     unsigned long indepVarCount = * _d_vector_res.rows() + {{2U|2UL|2ULL}};
+// CHECK-NEXT:     clad::array<double> _d_vector_i = clad::one_hot_vector(indepVarCount, {{0U|0UL|0ULL}});
+// CHECK-NEXT:     clad::array<double> _d_vector_j = clad::one_hot_vector(indepVarCount, {{1U|1UL|1ULL}});
+// CHECK-NEXT:     * _d_vector_res = clad::identity_matrix(* _d_vector_res.rows(), indepVarCount, {{2U|2UL|2ULL}});
+// CHECK-NEXT:     * _d_vector_res[0] = _d_vector_i * i + i * _d_vector_i;
 // CHECK-NEXT:     res[0] = i * i;
-// CHECK-NEXT:     double _t1 = res[1];
+// CHECK-NEXT:     * _d_vector_res[1] = _d_vector_j * j + j * _d_vector_j;
 // CHECK-NEXT:     res[1] = j * j;
-// CHECK-NEXT:     double _t2 = res[2];
+// CHECK-NEXT:     * _d_vector_res[2] = _d_vector_i * j + i * _d_vector_j;
 // CHECK-NEXT:     res[2] = i * j;
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{4U|4UL|4ULL}}] += 1 * j;
-// CHECK-NEXT:             jacobianMatrix[{{5U|5UL|5ULL}}] += i * 1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         res[2] = _t2;
-// CHECK-NEXT:     }
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += 1 * j;
-// CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += j * 1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         res[1] = _t1;
-// CHECK-NEXT:     }
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += 1 * i;
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += i * 1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         res[0] = _t0;
-// CHECK-NEXT:     }
 // CHECK-NEXT: }
 
 
@@ -55,41 +37,22 @@ constexpr void f_1(double x, double y, double z, double output[]) {
   output[2] = z * x * 10 - y * z;
 }
 
-// CHECK: constexpr void f_1_jac(double x, double y, double z, double output[], double *jacobianMatrix) {
-// CHECK-NEXT:     double _t0 = output[0];
-// CHECK-NEXT:     output[0] = x * x * x; 
-// CHECK-NEXT:     double _t1 = output[1];
-// CHECK-NEXT:     output[1] = x * y * x + y * x * x; 
-// CHECK-NEXT:     double _t2 = output[2];
-// CHECK-NEXT:     output[2] = z * x * 10 - y * z; 
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{8U|8UL|8ULL}}] += 1 * 10 * x; 
-// CHECK-NEXT:             jacobianMatrix[{{6U|6UL|6ULL}}] += z * 1 * 10; 
-// CHECK-NEXT:             jacobianMatrix[{{7U|7UL|7ULL}}] += -1 * z;     
-// CHECK-NEXT:             jacobianMatrix[{{8U|8UL|8ULL}}] += y * -1;     
-// CHECK-NEXT:         }
-// CHECK-NEXT:         output[2] = _t2; 
-// CHECK-NEXT:     }
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += 1 * x * y; 
-// CHECK-NEXT:             jacobianMatrix[{{4U|4UL|4ULL}}] += x * 1 * x;
-// CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += x * y * 1;
-// CHECK-NEXT:             jacobianMatrix[{{4U|4UL|4ULL}}] += 1 * x * x;
-// CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += y * 1 * x; 
-// CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += y * x * 1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         output[1] = _t1; 
-// CHECK-NEXT:     }
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += 1 * x * x; 
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += x * 1 * x; 
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += x * x * 1; 
-// CHECK-NEXT:         }
-// CHECK-NEXT:         output[0] = _t0;
-// CHECK-NEXT:     }
+// CHECK: constexpr void f_1_jac(double x, double y, double z, double output[], clad::array_ref<clad::matrix<double> > _d_vector_output) {
+// CHECK-NEXT:     unsigned long indepVarCount = * _d_vector_output.rows() + {{3U|3UL|3ULL}};
+// CHECK-NEXT:     clad::array<double> _d_vector_x = clad::one_hot_vector(indepVarCount, {{0U|0UL|0ULL}});
+// CHECK-NEXT:     clad::array<double> _d_vector_y = clad::one_hot_vector(indepVarCount, {{1U|1UL|1ULL}});
+// CHECK-NEXT:     clad::array<double> _d_vector_z = clad::one_hot_vector(indepVarCount, {{2U|2UL|2ULL}});
+// CHECK-NEXT:     * _d_vector_output = clad::identity_matrix(* _d_vector_output.rows(), indepVarCount, {{3U|3UL|3ULL}});
+// CHECK-NEXT:     double _t0 = x * x;
+// CHECK-NEXT:     * _d_vector_output[0] = (_d_vector_x * x + x * _d_vector_x) * x + _t0 * _d_vector_x;
+// CHECK-NEXT:     output[0] = _t0 * x;
+// CHECK-NEXT:     double _t1 = x * y;
+// CHECK-NEXT:     double _t2 = y * x;
+// CHECK-NEXT:     * _d_vector_output[1] = (_d_vector_x * y + x * _d_vector_y) * x + _t1 * _d_vector_x + (_d_vector_y * x + y * _d_vector_x) * x + _t2 * _d_vector_x;
+// CHECK-NEXT:     output[1] = _t1 * x + _t2 * x;
+// CHECK-NEXT:     double _t3 = z * x;
+// CHECK-NEXT:     * _d_vector_output[2] = (_d_vector_z * x + z * _d_vector_x) * 10 + _t3 * 0 - (_d_vector_y * z + y * _d_vector_z);
+// CHECK-NEXT:     output[2] = _t3 * 10 - y * z;
 // CHECK-NEXT: }
 
 int main() {
@@ -97,7 +60,7 @@ int main() {
     INIT_JACOBIAN(fn_mul);
     INIT_JACOBIAN(f_1);
 
-    TEST_JACOBIAN(fn_mul, 2, 6, 3, 1, result, jacobianou); // CHECK-EXEC: {6.00, 0.00, 0.00, 2.00, 1.00, 3.00}
-    TEST_JACOBIAN(f_1, 3, 9, 4, 5, 6, result1, jacobianou1); // CHECK-EXEC: {48.00, 0.00, 0.00, 80.00, 32.00, 0.00, 60.00, -6.00, 35.00}
+    TEST_JACOBIAN(fn_mul, 2, 6, 3, 1, result, &jacobian); // CHECK-EXEC: {6.00, 0.00, 0.00, 2.00, 1.00, 3.00}
+    TEST_JACOBIAN(f_1, 3, 9, 4, 5, 6, result1, &jacobian); // CHECK-EXEC: {48.00, 0.00, 0.00, 80.00, 32.00, 0.00, 60.00, -6.00, 35.00}
 }
 
diff --git a/test/Jacobian/testUtility.C b/test/Jacobian/testUtility.C
index 7b38eb7ec..5ec505394 100644
--- a/test/Jacobian/testUtility.C
+++ b/test/Jacobian/testUtility.C
@@ -8,9 +8,8 @@
 #include "../TestUtils.h"
 
   double output[3] = {0};
-  double jacobian[6] = {0};
   double output1[3] = {0};
-  double jacobian1[9] = {0};
+  clad::matrix<double> jacobian(3, 3);
 
 void fn_mul(double i, double j, double *res) {
    res[0] = i*i;
@@ -18,34 +17,17 @@ void fn_mul(double i, double j, double *res) {
    res[2] = i*j;
 }
 
-// CHECK: void fn_mul_jac(double i, double j, double *res, double *jacobianMatrix) {
-// CHECK-NEXT:     double _t0 = res[0];
+// CHECK: void fn_mul_jac(double i, double j, double *res, clad::array_ref<clad::matrix<double> > _d_vector_res) {
+// CHECK-NEXT:     unsigned long indepVarCount = * _d_vector_res.rows() + {{2U|2UL|2ULL}};
+// CHECK-NEXT:     clad::array<double> _d_vector_i = clad::one_hot_vector(indepVarCount, {{0U|0UL|0ULL}});
+// CHECK-NEXT:     clad::array<double> _d_vector_j = clad::one_hot_vector(indepVarCount, {{1U|1UL|1ULL}});
+// CHECK-NEXT:     * _d_vector_res = clad::identity_matrix(* _d_vector_res.rows(), indepVarCount, {{2U|2UL|2ULL}});
+// CHECK-NEXT:     * _d_vector_res[0] = _d_vector_i * i + i * _d_vector_i;
 // CHECK-NEXT:     res[0] = i * i;
-// CHECK-NEXT:     double _t1 = res[1];
+// CHECK-NEXT:     * _d_vector_res[1] = _d_vector_j * j + j * _d_vector_j;
 // CHECK-NEXT:     res[1] = j * j;
-// CHECK-NEXT:     double _t2 = res[2];
+// CHECK-NEXT:     * _d_vector_res[2] = _d_vector_i * j + i * _d_vector_j;
 // CHECK-NEXT:     res[2] = i * j;
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{4U|4UL|4ULL}}] += 1 * j;
-// CHECK-NEXT:             jacobianMatrix[{{5U|5UL|5ULL}}] += i * 1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         res[2] = _t2;
-// CHECK-NEXT:     }
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += 1 * j;
-// CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += j * 1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         res[1] = _t1;
-// CHECK-NEXT:     }
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += 1 * i;
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += i * 1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         res[0] = _t0;
-// CHECK-NEXT:     }
 // CHECK-NEXT: }
 
 
@@ -56,41 +38,22 @@ void f_1(double x, double y, double z, double output[]) {
   output[2] = z * x * 10 - y * z;
 }
 
-// CHECK: void f_1_jac(double x, double y, double z, double output[], double *jacobianMatrix) {
-// CHECK-NEXT:     double _t0 = output[0];
-// CHECK-NEXT:     output[0] = x * x * x;
-// CHECK-NEXT:     double _t1 = output[1];
-// CHECK-NEXT:     output[1] = x * y * x + y * x * x;
-// CHECK-NEXT:     double _t2 = output[2];
-// CHECK-NEXT:     output[2] = z * x * 10 - y * z;
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{8U|8UL|8ULL}}] += 1 * 10 * x;
-// CHECK-NEXT:             jacobianMatrix[{{6U|6UL|6ULL}}] += z * 1 * 10;
-// CHECK-NEXT:             jacobianMatrix[{{7U|7UL|7ULL}}] += -1 * z;
-// CHECK-NEXT:             jacobianMatrix[{{8U|8UL|8ULL}}] += y * -1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         output[2] = _t2;
-// CHECK-NEXT:     }
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += 1 * x * y;
-// CHECK-NEXT:             jacobianMatrix[{{4U|4UL|4ULL}}] += x * 1 * x;
-// CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3ULL}}] += x * y * 1;
-// CHECK-NEXT:             jacobianMatrix[{{4U|4UL|4ULL}}] += 1 * x * x;
-// CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3OLL}}] += y * 1 * x;
-// CHECK-NEXT:             jacobianMatrix[{{3U|3UL|3OLL}}] += y * x * 1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         output[1] = _t1;
-// CHECK-NEXT:     }
-// CHECK-NEXT:     {
-// CHECK-NEXT:         {
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += 1 * x * x;
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += x * 1 * x;
-// CHECK-NEXT:             jacobianMatrix[{{0U|0UL|0ULL}}] += x * x * 1;
-// CHECK-NEXT:         }
-// CHECK-NEXT:         output[0] = _t0;
-// CHECK-NEXT:     }
+// CHECK: void f_1_jac(double x, double y, double z, double output[], clad::array_ref<clad::matrix<double> > _d_vector_output) {
+// CHECK-NEXT:     unsigned long indepVarCount = * _d_vector_output.rows() + {{3U|3UL|3ULL}};
+// CHECK-NEXT:     clad::array<double> _d_vector_x = clad::one_hot_vector(indepVarCount, {{0U|0UL|0ULL}});
+// CHECK-NEXT:     clad::array<double> _d_vector_y = clad::one_hot_vector(indepVarCount, {{1U|1UL|1ULL}});
+// CHECK-NEXT:     clad::array<double> _d_vector_z = clad::one_hot_vector(indepVarCount, {{2U|2UL|2ULL}});
+// CHECK-NEXT:     * _d_vector_output = clad::identity_matrix(* _d_vector_output.rows(), indepVarCount, {{3U|3UL|3ULL}});
+// CHECK-NEXT:     double _t0 = x * x;
+// CHECK-NEXT:     * _d_vector_output[0] = (_d_vector_x * x + x * _d_vector_x) * x + _t0 * _d_vector_x;
+// CHECK-NEXT:     output[0] = _t0 * x;
+// CHECK-NEXT:     double _t1 = x * y;
+// CHECK-NEXT:     double _t2 = y * x;
+// CHECK-NEXT:     * _d_vector_output[1] = (_d_vector_x * y + x * _d_vector_y) * x + _t1 * _d_vector_x + (_d_vector_y * x + y * _d_vector_x) * x + _t2 * _d_vector_x;
+// CHECK-NEXT:     output[1] = _t1 * x + _t2 * x;
+// CHECK-NEXT:     double _t3 = z * x;
+// CHECK-NEXT:     * _d_vector_output[2] = (_d_vector_z * x + z * _d_vector_x) * 10 + _t3 * 0 - (_d_vector_y * z + y * _d_vector_z);
+// CHECK-NEXT:     output[2] = _t3 * 10 - y * z;
 // CHECK-NEXT: }
 
 
@@ -98,7 +61,7 @@ int main(){
     INIT_JACOBIAN(fn_mul);
     INIT_JACOBIAN(f_1);
 
-    TEST_JACOBIAN(fn_mul, 2, 6, 3, 1, output, jacobian); // CHECK-EXEC: {6.00, 0.00, 0.00, 2.00, 1.00, 3.00}
-    TEST_JACOBIAN(f_1, 3, 9, 4, 5, 6, output1, jacobian1); // CHECK-EXEC: {48.00, 0.00, 0.00, 80.00, 32.00, 0.00, 60.00, -6.00, 35.00}
+    TEST_JACOBIAN(fn_mul, 2, 6, 3, 1, output, &jacobian); // CHECK-EXEC: {6.00, 0.00, 0.00, 2.00, 1.00, 3.00}
+    TEST_JACOBIAN(f_1, 3, 9, 4, 5, 6, output1, &jacobian); // CHECK-EXEC: {48.00, 0.00, 0.00, 80.00, 32.00, 0.00, 60.00, -6.00, 35.00}
 
 }
diff --git a/test/TestUtils.h b/test/TestUtils.h
index 6d76c59f6..13d78df2c 100644
--- a/test/TestUtils.h
+++ b/test/TestUtils.h
@@ -116,16 +116,25 @@ void run_gradient_impl(CF cf, index_pack<S...> s, Args&&... args) {
   display(std::get<S>(t)...);
 }
 
-template <class CF, std::size_t... S, class... Args>
-void run_jacobian_impl(CF cf, std::size_t size, index_pack<S...> s,
-                       Args&&... args) {
+template <class CF, class... Args>
+void run_jacobian_impl(CF cf, std::size_t size, Args&&... args) {
   std::tuple<Args...> t = {args...};
-  reset(std::get<S>(t)...);
   cf.execute(args...);
   typedef
       typename std::remove_reference<decltype(*std::get<sizeof...(args) - 1>(
           t))>::type arrElemType;
-  displayarr<arrElemType>(std::get<sizeof...(args) - 1>(t), size);
+  auto& res = *std::get<sizeof...(args) - 1>(t);
+  unsigned numOfParameters = sizeof...(args) - 2;
+  unsigned numOfOutputs = res.rows();
+  printf("{");
+  for (unsigned i = 0; i < numOfOutputs; ++i) {
+    for (unsigned j = 0; j < numOfParameters; ++j) {
+      printf("%.2f", res[i][j]);
+      if (i != numOfOutputs - 1 || j != numOfParameters - 1)
+        printf(", ");
+    }
+  }
+  printf("}\n");
 }
 
 template <class CF, std::size_t... S, class... Args>
@@ -150,11 +159,7 @@ void run_gradient(CF cf, Args&&... args) {
 template <std::size_t NumOfDerivativeArgs, std::size_t size, class CF,
           class... Args>
 void run_jacobian(CF cf, Args&&... args) {
-  using DerivativeArgRange =
-      typename GenerateRange<sizeof...(Args) - NumOfDerivativeArgs,
-                             sizeof...(Args) - 1>::type;
-  run_jacobian_impl(cf, size, DerivativeArgRange(),
-                    std::forward<Args>(args)...);
+  run_jacobian_impl(cf, size, std::forward<Args>(args)...);
 }
 
 template <std::size_t NumOfDerivativeArgs, class CF, class... Args>