diff --git a/CMakeLists.txt b/CMakeLists.txt
index 29b8cf90..7428b046 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -103,6 +103,8 @@ endif()
 
 if(${GRIDKIT_ENABLE_ENZYME})
   include(FindEnzyme)
+  # todo Add a centralized configuration file
+  add_definitions(-DGRIDKIT_ENABLE_ENZYME)
 endif()
 
 # Macro that adds libraries
diff --git a/cmake/EnzymeAddLibrary.cmake b/cmake/EnzymeAddLibrary.cmake
new file mode 100644
index 00000000..9d7832ee
--- /dev/null
+++ b/cmake/EnzymeAddLibrary.cmake
@@ -0,0 +1,123 @@
+# 
+#[[
+
+Macro to manually compile with Enzyme
+
+Author(s):
+- Asher Mancinelli <ashermancinelli@gmail.com>
+- Nicholson Koukpaizan <koukpaizannk@ornl.gov>
+
+]]
+
+macro(enzyme_build_object)
+  set(options)
+  set(oneValueArgs NAME)
+  set(multiValueArgs SOURCES LINK_LIBRARIES INCLUDE_DIRECTORIES)
+  cmake_parse_arguments(enzyme_build_object "${options}" "${oneValueArgs}"
+    "${multiValueArgs}" ${ARGN})
+
+  set(PHASE2 "${CMAKE_CURRENT_BINARY_DIR}/${enzyme_build_object_NAME}.bc")
+  set(PHASE3 "${CMAKE_CURRENT_BINARY_DIR}/${enzyme_build_object_NAME}_enzyme.ll")
+  set(PHASE4 "${CMAKE_CURRENT_BINARY_DIR}/${enzyme_build_object_NAME}_opt.ll")
+  set(PHASE5 "${CMAKE_CURRENT_BINARY_DIR}/${enzyme_build_object_NAME}")
+
+  set(OBJS "")
+  set(includes "${enzyme_build_object_INCLUDE_DIRECTORIES}")
+
+  foreach(lib ${enzyme_build_object_LINK_LIBRARIES})
+    get_target_property(include ${lib} INCLUDE_DIRECTORIES)
+    set(includes "${includes}" ${include})
+
+    get_target_property(libsource ${lib} SOURCES)
+    string(FIND "${libsource}" "TARGET" found)
+    if(NOT(${found} EQUAL -1))
+      list(APPEND LINKER_FLAGS "-Wl,${libsource}")
+    endif()
+  endforeach()
+
+  foreach(dir ${includes})
+    if(EXISTS ${dir})
+      list(APPEND INCLUDE_COMPILER_LIST "-I${dir}")
+    endif()
+  endforeach()
+
+  foreach(SRC ${enzyme_build_object_SOURCES})
+    set(PHASE0 "${CMAKE_CURRENT_SOURCE_DIR}/${SRC}")
+    set(PHASE1 "${CMAKE_CURRENT_BINARY_DIR}/${enzyme_build_object_NAME}_${SRC}_compile.o")
+    add_custom_command(
+      DEPENDS ${PHASE0} 
+      OUTPUT ${PHASE1}
+      COMMAND ${CMAKE_CXX_COMPILER} -flto -c ${PHASE0} ${INCLUDE_COMPILER_LIST} -O2 -fno-vectorize -ffast-math -fno-unroll-loops -fpass-plugin=${ENZYME_CLANG_PLUGIN_LIBRARY} -Xclang -load -Xclang ${ENZYME_CLANG_PLUGIN_LIBRARY} -mllvm -enable-load-pre=0 -mllvm -enzyme-auto-sparsity=1 -o ${PHASE1}
+      COMMENT "Compiling ${SRC} to object file for target ${enzyme_build_object_NAME}"
+      )
+    set(OBJS "${OBJS} ${PHASE1}")
+  endforeach()
+
+  cmake_language(EVAL CODE "
+  add_custom_command(
+    DEPENDS ${OBJS}
+    OUTPUT ${PHASE2}
+    COMMAND ${GRIDKIT_LLVM_LINK} ${OBJS} -o ${PHASE2}
+    COMMENT \"Linking object files to LLVM bytecode for target ${enzyme_build_object_NAME}\"
+    )
+  ")
+
+  add_custom_command(
+    DEPENDS ${PHASE2}
+    OUTPUT ${PHASE3}
+    COMMAND ${GRIDKIT_OPT} ${PHASE2} -load-pass-plugin=${ENZYME_LLVM_PLUGIN_LIBRARY} -passes=enzyme -o ${PHASE3} -S
+    COMMENT "Running Enzyme opt pass on target ${enzyme_build_object_NAME}"
+    )
+
+  add_custom_command(
+    DEPENDS ${PHASE3}
+    OUTPUT ${PHASE4}
+    COMMAND ${GRIDKIT_OPT} ${PHASE3} -O2 -o ${PHASE4} -S
+    COMMENT "Running remaining opt passes on target ${enzyme_build_object_NAME}"
+    )
+
+  add_custom_command(
+    DEPENDS ${PHASE4}
+    OUTPUT ${PHASE5}
+    COMMAND ${CMAKE_CXX_COMPILER} -c ${PHASE4} -o ${PHASE5}
+    COMMENT "Generating optimized object file for target ${enzyme_build_object_NAME}"
+    )
+endmacro()
+
+macro(enzyme_add_executable)
+  set(options)
+  set(oneValueArgs NAME)
+  set(multiValueArgs SOURCES LINK_LIBRARIES INCLUDE_DIRECTORIES)
+  cmake_parse_arguments(enzyme_add_executable "${options}" "${oneValueArgs}"
+    "${multiValueArgs}" ${ARGN})
+
+  enzyme_build_object(
+    NAME "${enzyme_add_executable_NAME}.o"
+    SOURCES ${enzyme_add_executable_SOURCES}
+    LINK_LIBRARIES ${enzyme_add_executable_LINK_LIBRARIES}
+    INCLUDE_DIRECTORIES ${enzyme_add_executable_INCLUDE_DIRECTORIES}
+  )
+
+  add_executable("${enzyme_add_executable_NAME}" "${enzyme_add_executable_NAME}.o")
+  set_target_properties("${enzyme_add_executable_NAME}" PROPERTIES LINKER_LANGUAGE CXX)
+  target_link_libraries("${enzyme_add_executable_NAME}" ${enzyme_add_executable_LINK_LIBRARIES})
+endmacro()
+
+macro(enzyme_add_library)
+  set(options)
+  set(oneValueArgs NAME)
+  set(multiValueArgs SOURCES LINK_LIBRARIES INCLUDE_DIRECTORIES)
+  cmake_parse_arguments(enzyme_add_library "${options}" "${oneValueArgs}"
+    "${multiValueArgs}" ${ARGN})
+
+  enzyme_build_object(
+    NAME "${enzyme_add_library_NAME}.o"
+    SOURCES ${enzyme_add_library_SOURCES}
+    LINK_LIBRARIES ${enzyme_add_library_LINK_LIBRARIES}
+    INCLUDE_DIRECTORIES ${enzyme_add_library_INCLUDE_DIRECTORIES}
+  )
+
+  add_library("${enzyme_add_library_NAME}" "${enzyme_add_library_NAME}.o")
+  set_target_properties("${enzyme_add_library_NAME}" PROPERTIES LINKER_LANGUAGE CXX)
+  target_link_libraries("${enzyme_add_library_NAME}" ${enzyme_add_library_LINK_LIBRARIES})
+endmacro()
diff --git a/cmake/FindEnzyme.cmake b/cmake/FindEnzyme.cmake
index 00a76d37..fbfaa440 100644
--- a/cmake/FindEnzyme.cmake
+++ b/cmake/FindEnzyme.cmake
@@ -59,82 +59,3 @@ find_program(GRIDKIT_OPT opt
              bin
              REQUIRED)
 message(STATUS "opt: ${GRIDKIT_OPT}")
-
-macro(enzyme_add_executable)
-  set(options)
-  set(oneValueArgs NAME)
-  set(multiValueArgs SOURCES LINK_LIBRARIES INCLUDE_DIRECTORIES)
-  cmake_parse_arguments(enzyme_add_executable "${options}" "${oneValueArgs}"
-    "${multiValueArgs}" ${ARGN})
-
-  set(PHASE2 "${CMAKE_CURRENT_BINARY_DIR}/${enzyme_add_executable_NAME}.bc")
-  set(PHASE3 "${CMAKE_CURRENT_BINARY_DIR}/${enzyme_add_executable_NAME}_enzyme.ll")
-  set(PHASE4 "${CMAKE_CURRENT_BINARY_DIR}/${enzyme_add_executable_NAME}_opt.ll")
-  set(PHASE5 "${CMAKE_CURRENT_BINARY_DIR}/${enzyme_add_executable_NAME}")
-
-  set(OBJS "")
-  set(includes "${enzyme_add_executable_INCLUDE_DIRECTORIES}")
-
-  foreach(lib ${enzyme_add_executable_LINK_LIBRARIES})
-    get_target_property(include ${lib} INCLUDE_DIRECTORIES)
-    set(includes "${includes}" ${include})
-
-    get_target_property(libsource ${lib} SOURCES)
-    string(FIND "${libsource}" "TARGET" found)
-    if(NOT(${found} EQUAL -1))
-      list(APPEND LINKER_FLAGS "-Wl,${libsource}")
-    endif()
-  endforeach()
-
-  foreach(dir ${includes})
-    if(EXISTS ${dir})
-      list(APPEND INCLUDE_COMPILER_LIST "-I${dir}")
-    endif()
-  endforeach()
-
-  foreach(SRC ${enzyme_add_executable_SOURCES})
-    set(PHASE0 "${CMAKE_CURRENT_SOURCE_DIR}/${SRC}")
-    set(PHASE1 "${CMAKE_CURRENT_BINARY_DIR}/${enzyme_add_executable_NAME}_${SRC}_compile.o")
-    add_custom_command(
-      DEPENDS ${PHASE0} 
-      OUTPUT ${PHASE1}
-      COMMAND ${CMAKE_CXX_COMPILER} -flto -c ${PHASE0} ${INCLUDE_COMPILER_LIST} -O2 -fno-vectorize -ffast-math -fno-unroll-loops -fpass-plugin=${ENZYME_CLANG_PLUGIN_LIBRARY} -Xclang -load -Xclang ${ENZYME_CLANG_PLUGIN_LIBRARY} -mllvm -enable-load-pre=0 -mllvm -enzyme-auto-sparsity=1 -o ${PHASE1}
-      COMMENT "Compiling ${SRC} to object file for target ${enzyme_add_executable_NAME}"
-      )
-    set(OBJS "${OBJS} ${PHASE1}")
-  endforeach()
-
-  cmake_language(EVAL CODE "
-  add_custom_command(
-    DEPENDS ${OBJS}
-    OUTPUT ${PHASE2}
-    COMMAND ${GRIDKIT_LLVM_LINK} ${OBJS} -o ${PHASE2}
-    COMMENT \"Linking object files to LLVM bytecode for target ${enzyme_add_executable_NAME}\"
-    )
-  ")
-
-  add_custom_command(
-    DEPENDS ${PHASE2}
-    OUTPUT ${PHASE3}
-    COMMAND ${GRIDKIT_OPT} ${PHASE2} -load-pass-plugin=${ENZYME_LLVM_PLUGIN_LIBRARY} -passes=enzyme -o ${PHASE3} -S
-    COMMENT "Running Enzyme opt pass on target ${enzyme_add_executable_NAME}"
-    )
-
-  add_custom_command(
-    DEPENDS ${PHASE3}
-    OUTPUT ${PHASE4}
-    COMMAND ${GRIDKIT_OPT} ${PHASE3} -O2 -o ${PHASE4} -S
-    COMMENT "Running remaining opt passes on target ${enzyme_add_executable_NAME}"
-    )
-
-  add_custom_command(
-    DEPENDS ${PHASE4} ${enzyme_add_executable_LINK_LIBRARIES}
-    OUTPUT ${PHASE5}
-    COMMAND ${CMAKE_CXX_COMPILER} ${LINKER_FLAGS} ${PHASE4} -o ${PHASE5}
-    )
-
-  add_custom_target(
-    "${enzyme_add_executable_NAME}_target" ALL
-    DEPENDS ${PHASE5}
-    )
-endmacro()
diff --git a/examples/Enzyme/Library/Scalar/CMakeLists.txt b/examples/Enzyme/Library/Scalar/CMakeLists.txt
index 6dff430a..21e1b9b9 100644
--- a/examples/Enzyme/Library/Scalar/CMakeLists.txt
+++ b/examples/Enzyme/Library/Scalar/CMakeLists.txt
@@ -1,6 +1,4 @@
-enzyme_add_executable(
-  NAME EnzymeLibScalarCheck 
-  SOURCES EnzymeScalar.cpp ScalarModel.cpp
-)
+add_executable(EnzymeLibScalarCheck EnzymeScalar.cpp ScalarModel.cpp)
+target_link_libraries(EnzymeLibScalarCheck ClangEnzymeFlags GRIDKIT::Utilities)
 
 add_test(NAME "EnzymeLibScalarCheck" COMMAND ${CMAKE_CURRENT_BINARY_DIR}/EnzymeLibScalarCheck)
diff --git a/examples/Enzyme/Library/Scalar/EnzymeScalar.cpp b/examples/Enzyme/Library/Scalar/EnzymeScalar.cpp
index a20f69e1..933e1bd7 100644
--- a/examples/Enzyme/Library/Scalar/EnzymeScalar.cpp
+++ b/examples/Enzyme/Library/Scalar/EnzymeScalar.cpp
@@ -1,7 +1,7 @@
 #include <iostream>
-#include <limits>
 
 #include "ScalarModel.hpp"
+#include <Utilities/Testing.hpp>
 
 /**
  * @brief Example that computes the derivative of a library function
@@ -23,7 +23,7 @@ int main()
   double dsq = scalar_model.getDerivativeValue();
 
   std::cout << "x = " << var << ", x^2 = " << sq << ", d(x^2)/dx = " << dsq << "\n";
-  if (std::abs(dsq - 2.0 * var) > std::numeric_limits<double>::epsilon())
+  if (!GridKit::Testing::isEqual(dsq, 2.0 * var))
   {
     fail++;
     std::cout << "Result incorrect\n";
diff --git a/examples/Enzyme/Library/Vector/CMakeLists.txt b/examples/Enzyme/Library/Vector/CMakeLists.txt
index eda952de..559a249b 100644
--- a/examples/Enzyme/Library/Vector/CMakeLists.txt
+++ b/examples/Enzyme/Library/Vector/CMakeLists.txt
@@ -1,7 +1,4 @@
-enzyme_add_executable(
-  NAME EnzymeLibVectorCheck 
-  SOURCES EnzymeVector.cpp VectorModel.cpp
-  LINK_LIBRARIES GRIDKIT::DenseMatrix
-)
+add_executable(EnzymeLibVectorCheck EnzymeVector.cpp VectorModel.cpp)
+target_link_libraries(EnzymeLibVectorCheck ClangEnzymeFlags GRIDKIT::DenseMatrix GRIDKIT::Utilities)
 
 add_test(NAME "EnzymeLibVectorCheck" COMMAND ${CMAKE_CURRENT_BINARY_DIR}/EnzymeLibVectorCheck)
diff --git a/examples/Enzyme/Library/Vector/EnzymeVector.cpp b/examples/Enzyme/Library/Vector/EnzymeVector.cpp
index 8a3e6925..4d9dd22f 100644
--- a/examples/Enzyme/Library/Vector/EnzymeVector.cpp
+++ b/examples/Enzyme/Library/Vector/EnzymeVector.cpp
@@ -1,7 +1,7 @@
 #include <iostream>
-#include <limits>
 
 #include "VectorModel.hpp"
+#include <Utilities/Testing.hpp>
 
 /**
  * @brief Example that computes the Jacobian of a vector-valued residual
@@ -20,12 +20,12 @@ inline double dsquare_ref_scalar(double x)
 DenseMatrix dsquare_ref(std::vector<double> x, std::vector<double> y)
 {
   DenseMatrix jac(x.size(), y.size());
-  for (int idy = 0; idy < y.size(); ++idy)
+  for (size_t idy = 0; idy < y.size(); ++idy)
   {
-    for (int idx = 0; idx < x.size(); ++idx)
+    for (size_t idx = 0; idx < x.size(); ++idx)
     {
-      if (idx == idy)
-        jac.setValue(idx, idy, dsquare_ref_scalar(x[idx]));
+      if (idy <= idx)
+        jac.setValue(idx, idy, dsquare_ref_scalar(x[idy]));
     }
   }
   return jac;
@@ -34,12 +34,12 @@ DenseMatrix dsquare_ref(std::vector<double> x, std::vector<double> y)
 int main()
 {
   // Size and variable declarations
-  constexpr int       n = 10;
+  constexpr size_t    n = 10;
   std::vector<double> var(n);
 
   // Random input values
-  srand(time(NULL));
-  for (int idx = 0; idx < var.size(); ++idx)
+  srand(static_cast<unsigned int>(time(NULL)));
+  for (size_t idx = 0; idx < var.size(); ++idx)
   {
     var[idx] = rand();
   }
@@ -59,11 +59,11 @@ int main()
   // Check
   int  fail    = 0;
   bool verbose = true;
-  for (int idy = 0; idy < res.size(); ++idy)
+  for (size_t idy = 0; idy < res.size(); ++idy)
   {
-    for (int idx = 0; idx < var.size(); ++idx)
+    for (size_t idx = 0; idx < var.size(); ++idx)
     {
-      if (std::abs(jac.getValue(idx, idy) - jac_ref.getValue(idx, idy)) > std::numeric_limits<double>::epsilon())
+      if (!GridKit::Testing::isEqual(jac.getValue(idx, idy), jac_ref.getValue(idx, idy)))
       {
         fail++;
         if (verbose)
diff --git a/examples/Enzyme/Library/Vector/VectorModel.cpp b/examples/Enzyme/Library/Vector/VectorModel.cpp
index cb4d89f5..bdd99842 100644
--- a/examples/Enzyme/Library/Vector/VectorModel.cpp
+++ b/examples/Enzyme/Library/Vector/VectorModel.cpp
@@ -4,7 +4,7 @@
 
 #include "EnzymeWrapper.hpp"
 
-VectorModel::VectorModel(int n)
+VectorModel::VectorModel(size_t n)
   : x_(n),
     f_(n),
     df_dx_(n, n)
@@ -18,15 +18,19 @@ inline double VectorModel::square_scalar(double x)
 
 void VectorModel::square(std::vector<double>& x, std::vector<double>& y)
 {
-  for (int idx = 0; idx < x.size(); ++idx)
+  for (size_t idx = 0; idx < x.size(); ++idx)
   {
-    y[idx] = this->square_scalar(x[idx]);
+    y[idx] = 0.0;
+    for (size_t idy = 0; idy <= idx; idy++)
+    {
+      y[idx] += this->square_scalar(x[idy]);
+    }
   }
 }
 
 void VectorModel::setVariable(std::vector<double> x)
 {
-  for (int idx = 0; idx < x.size(); ++idx)
+  for (size_t idx = 0; idx < x.size(); ++idx)
   {
     x_[idx] = x[idx];
   }
@@ -39,13 +43,13 @@ void VectorModel::evalResidual()
 
 void VectorModel::evalJacobian()
 {
-  const int           n = x_.size();
+  const size_t        n = x_.size();
   std::vector<double> v(n);
   VectorModel         d_vector_model(n);
-  for (int idy = 0; idy < n; ++idy)
+  for (size_t idy = 0; idy < n; ++idy)
   {
     // Elementary vector for Jacobian-vector product
-    for (int idx = 0; idx < n; ++idx)
+    for (size_t idx = 0; idx < n; ++idx)
     {
       v[idx] = 0.0;
     }
@@ -60,7 +64,7 @@ void VectorModel::evalJacobian()
         &d_vector_model);
 
     // Store result
-    for (int idx = 0; idx < n; ++idx)
+    for (size_t idx = 0; idx < n; ++idx)
     {
       df_dx_.setValue(idx, idy, d_res[idx]);
     }
diff --git a/examples/Enzyme/Library/Vector/VectorModel.hpp b/examples/Enzyme/Library/Vector/VectorModel.hpp
index da66f696..fd85f716 100644
--- a/examples/Enzyme/Library/Vector/VectorModel.hpp
+++ b/examples/Enzyme/Library/Vector/VectorModel.hpp
@@ -17,7 +17,7 @@ class VectorModel
   void                square(std::vector<double>&, std::vector<double>&);
 
 public:
-  VectorModel(int);
+  VectorModel(size_t);
   void                 setVariable(std::vector<double>);
   void                 evalResidual();
   void                 evalJacobian();
diff --git a/examples/Enzyme/PowerElectronics/CMakeLists.txt b/examples/Enzyme/PowerElectronics/CMakeLists.txt
index a8e74e88..a24c2326 100644
--- a/examples/Enzyme/PowerElectronics/CMakeLists.txt
+++ b/examples/Enzyme/PowerElectronics/CMakeLists.txt
@@ -1,7 +1,5 @@
-enzyme_add_executable(
-  NAME EnzymePowerElectronicsCheck 
-  SOURCES main.cpp
-  LINK_LIBRARIES GRIDKIT::DenseMatrix GRIDKIT::power_elec_disgen
-)
+add_executable(EnzymePowerElectronicsCheck main.cpp)
+target_compile_options(EnzymePowerElectronicsCheck PUBLIC -fno-vectorize -ffast-math -fno-unroll-loops)
+target_link_libraries(EnzymePowerElectronicsCheck ClangEnzymeFlags GRIDKIT::DenseMatrix GRIDKIT::power_elec_disgen GRIDKIT::Utilities)
 
 add_test(NAME "EnzymePowerElectronicsCheck" COMMAND ${CMAKE_CURRENT_BINARY_DIR}/EnzymePowerElectronicsCheck)
diff --git a/examples/Enzyme/PowerElectronics/main.cpp b/examples/Enzyme/PowerElectronics/main.cpp
index 5698901d..35de5af9 100644
--- a/examples/Enzyme/PowerElectronics/main.cpp
+++ b/examples/Enzyme/PowerElectronics/main.cpp
@@ -1,9 +1,9 @@
 #include <iostream>
-#include <limits>
 
 #include <LinearAlgebra/DenseMatrix/DenseMatrix.hpp>
 #include <Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp>
 #include <Model/PowerElectronics/SystemModelPowerElectronics.hpp>
+#include <Utilities/Testing.hpp>
 
 /**
  * @brief Standalone example that computes the Jacobian associated with the
@@ -97,15 +97,15 @@ void evaluateResidual(std::vector<double> y_, std::vector<double> f_)
 template <typename T>
 void EnzymeModelJacobian(T* model, DenseMatrix& jac)
 {
-  int                 N = model->size();
+  size_t              N = model->size();
   std::vector<double> y(N);
   std::vector<double> v(N);
   std::vector<double> res(N);
   std::vector<double> d_res(N);
-  for (int idy = 0; idy < N; ++idy)
+  for (size_t idy = 0; idy < N; ++idy)
   {
     // Elementary vector for Jacobian-vector product
-    for (int idx = 0; idx < N; ++idx)
+    for (size_t idx = 0; idx < N; ++idx)
     {
       y[idx]   = (model->y())[idx];
       res[idx] = (model->getResidual())[idx];
@@ -123,7 +123,7 @@ void EnzymeModelJacobian(T* model, DenseMatrix& jac)
                      &d_res);
 
     // Store result
-    for (int idx = 0; idx < N; ++idx)
+    for (size_t idx = 0; idx < N; ++idx)
     {
       jac.setValue(idx, idy, d_res[idx]);
     }
@@ -171,16 +171,18 @@ int main()
   // Check
   int  fail    = 0;
   bool verbose = true;
-  for (int idy = 0; idy < dg->size(); ++idy)
+  for (size_t idy = 0; idy < dg->size(); ++idy)
   {
-    for (int idx = 0; idx < dg->size(); ++idx)
+    for (size_t idx = 0; idx < dg->size(); ++idx)
     {
-      if (std::abs(jac_autodiff.getValue(idx, idy) - jac_ref_dense.getValue(idx, idy)) > std::numeric_limits<double>::epsilon())
+      double jac_value     = jac_autodiff.getValue(idx, idy);
+      double jac_ref_value = jac_ref_dense.getValue(idx, idy);
+      if (!GridKit::Testing::isEqual(jac_value, jac_ref_value))
       {
         fail++;
         if (verbose)
         {
-          std::cout << "Result incorrect at line = " << idy << ", column = " << idx << "\n";
+          std::cout << "Result incorrect at line = " << idy << ", column = " << idx << ", obtained = " << jac_value << ", reference = " << jac_ref_value << ", difference = " << std::abs(jac_value - jac_ref_value) << "\n";
         }
       }
     }
diff --git a/examples/Enzyme/Standalone/CMakeLists.txt b/examples/Enzyme/Standalone/CMakeLists.txt
index d80ba4f5..e62a2ea1 100644
--- a/examples/Enzyme/Standalone/CMakeLists.txt
+++ b/examples/Enzyme/Standalone/CMakeLists.txt
@@ -1,19 +1,12 @@
-enzyme_add_executable(
-  NAME EnzymeStandaloneScalarCheck
-  SOURCES EnzymeScalar.cpp
-)
+add_executable(EnzymeStandaloneScalarCheck EnzymeScalar.cpp)
+target_link_libraries(EnzymeStandaloneScalarCheck ClangEnzymeFlags GRIDKIT::Utilities)
 
-enzyme_add_executable(
-  NAME EnzymeStandaloneVectorCheck
-  SOURCES EnzymeVector.cpp
-  LINK_LIBRARIES GRIDKIT::DenseMatrix
-)
+add_executable(EnzymeStandaloneVectorCheck EnzymeVector.cpp)
+target_link_libraries(EnzymeStandaloneVectorCheck ClangEnzymeFlags GRIDKIT::DenseMatrix GRIDKIT::Utilities)
 
-enzyme_add_executable(
-  NAME EnzymeStandaloneSparseCheck
-  SOURCES EnzymeSparse.cpp
-  LINK_LIBRARIES GRIDKIT::SparseMatrix
-)
+add_executable(EnzymeStandaloneSparseCheck EnzymeSparse.cpp)
+target_compile_options(EnzymeStandaloneSparseCheck PUBLIC -mllvm -enzyme-auto-sparsity=1)
+target_link_libraries(EnzymeStandaloneSparseCheck ClangEnzymeFlags GRIDKIT::SparseMatrix GRIDKIT::Utilities)
 
 add_test(NAME "EnzymeStandaloneScalarCheck" COMMAND ${CMAKE_CURRENT_BINARY_DIR}/EnzymeStandaloneScalarCheck)
 add_test(NAME "EnzymeStandaloneVectorCheck" COMMAND ${CMAKE_CURRENT_BINARY_DIR}/EnzymeStandaloneVectorCheck)
diff --git a/examples/Enzyme/Standalone/EnzymeScalar.cpp b/examples/Enzyme/Standalone/EnzymeScalar.cpp
index d799e9b0..01f1ce8d 100644
--- a/examples/Enzyme/Standalone/EnzymeScalar.cpp
+++ b/examples/Enzyme/Standalone/EnzymeScalar.cpp
@@ -1,5 +1,6 @@
 #include <iostream>
-#include <limits>
+
+#include <Utilities/Testing.hpp>
 
 /**
  * @brief Standalone example that computes the derivative of a scalar function
@@ -27,7 +28,7 @@ int main()
   double sq   = square(var);
   double dsq  = dsquare(var);
   std::cout << "x = " << var << ", x^2 = " << sq << ", d(x^2)/dx = " << dsq << "\n";
-  if (std::abs(dsq - 2.0 * var) > std::numeric_limits<double>::epsilon())
+  if (!GridKit::Testing::isEqual(dsq, 2.0 * var))
   {
     fail++;
     std::cout << "Result incorrect\n";
diff --git a/examples/Enzyme/Standalone/EnzymeSparse.cpp b/examples/Enzyme/Standalone/EnzymeSparse.cpp
index f8e0bf0a..22311318 100644
--- a/examples/Enzyme/Standalone/EnzymeSparse.cpp
+++ b/examples/Enzyme/Standalone/EnzymeSparse.cpp
@@ -1,10 +1,10 @@
 #include <assert.h>
 #include <cmath>
-#include <limits>
 #include <stdio.h>
 #include <vector>
 
 #include <LinearAlgebra/SparseMatrix/COO_Matrix.hpp>
+#include <Utilities/Testing.hpp>
 
 /**
  * @brief Standalone example that computes the sparse Jacobian of a vector-valued function
@@ -26,15 +26,15 @@ template <typename T, typename... Tys>
 extern T __enzyme_todense(Tys...) noexcept;
 
 /// Sparse storage for Enzyme
-template <typename T>
+template <typename ScalarT>
 struct Triple
 {
-  size_t row;
-  size_t col;
-  T      val;
+  size_t  row;
+  size_t  col;
+  ScalarT val;
   Triple(Triple&&) = default;
 
-  Triple(size_t row, size_t col, T val)
+  Triple(size_t row, size_t col, ScalarT val)
     : row(row),
       col(col),
       val(val)
@@ -42,73 +42,77 @@ struct Triple
   }
 };
 
-__attribute__((enzyme_sparse_accumulate)) static void inner_storeflt(int64_t row, int64_t col, float val, std::vector<Triple<float>>& triplets)
+[[maybe_unused]] __attribute__((enzyme_sparse_accumulate)) static void inner_storeflt(size_t row, size_t col, float val, std::vector<Triple<float>>& triplets)
 {
   triplets.emplace_back(row, col, val);
 }
 
-__attribute__((enzyme_sparse_accumulate)) static void inner_storedbl(int64_t row, int64_t col, double val, std::vector<Triple<double>>& triplets)
+[[maybe_unused]] __attribute__((enzyme_sparse_accumulate)) static void inner_storedbl(size_t row, size_t col, double val, std::vector<Triple<double>>& triplets)
 {
   triplets.emplace_back(row, col, val);
 }
 
-template <typename T>
-__attribute__((always_inline)) static void sparse_store(T val, int64_t idx, size_t i, std::vector<Triple<T>>& triplets)
+template <typename ScalarT>
+__attribute__((always_inline)) static void sparse_store(ScalarT val, size_t idx, size_t i, std::vector<Triple<ScalarT>>& triplets)
 {
   if (val == 0.0)
     return;
-  idx /= sizeof(T);
-  if constexpr (sizeof(T) == 4)
-    inner_storeflt(i, idx, val, triplets);
+  idx /= sizeof(ScalarT);
+  if constexpr (sizeof(ScalarT) == 4)
+    inner_storeflt(idx, i, val, triplets);
   else
-    inner_storedbl(i, idx, val, triplets);
+    inner_storedbl(idx, i, val, triplets);
 }
 
-template <typename T>
-__attribute__((always_inline)) static T sparse_load(int64_t idx, size_t i, std::vector<Triple<T>>& triplets)
+template <typename ScalarT>
+__attribute__((always_inline)) static ScalarT sparse_load(size_t, size_t, std::vector<Triple<ScalarT>>&)
 {
   return 0.0;
 }
 
-template <typename T>
-__attribute__((always_inline)) static void ident_store(T, int64_t idx, size_t i)
+template <typename ScalarT>
+__attribute__((always_inline)) static void ident_store(ScalarT, size_t, size_t)
 {
   assert(0 && "should never load");
 }
 
-template <typename T>
-__attribute__((always_inline)) static T ident_load(int64_t idx, size_t i)
+template <typename ScalarT>
+__attribute__((always_inline)) static ScalarT ident_load(size_t idx, size_t i)
 {
-  idx /= sizeof(T);
-  return (T) (idx == i);
+  idx /= sizeof(ScalarT);
+  return (ScalarT) (idx == i);
 }
 
 /// Vector-valued function to differentiate
-template <typename T>
-__attribute__((always_inline)) static void f(size_t N, T* input, T* output)
+template <typename ScalarT>
+__attribute__((always_inline)) static void f(size_t N, ScalarT* input, ScalarT* output)
 {
-  for (size_t i = 0; i < N; i++)
+  for (size_t idx = 0; idx < N; ++idx)
   {
-    output[i] = input[i] * input[i];
+    output[idx] = 0.0;
+    for (size_t idy = 0; idy <= idx; idy++)
+    {
+      output[idx] += input[idy] * input[idy];
+    }
   }
 }
 
 /// Reference Jacobian
-template <typename T>
-void jac_f_ref(std::vector<T> x, std::vector<T> y, SparseMatrix& jac)
+template <typename ScalarT>
+void jac_f_ref(std::vector<ScalarT> x, std::vector<ScalarT> y, SparseMatrix& jac)
 {
-  std::vector<size_t> ctemp{};
-  std::vector<size_t> rtemp{};
-  std::vector<T>      valtemp{};
-  for (int idy = 0; idy < y.size(); ++idy)
+  std::vector<size_t>  ctemp{};
+  std::vector<size_t>  rtemp{};
+  std::vector<ScalarT> valtemp{};
+  for (size_t idy = 0; idy < y.size(); ++idy)
   {
-    for (int idx = 0; idx < x.size(); ++idx)
+    for (size_t idx = 0; idx < x.size(); ++idx)
     {
-      if (idx == idy)
+      if (idy <= idx)
       {
         rtemp.push_back(idx);
         ctemp.push_back(idy);
-        valtemp.push_back(2.0 * x[idx]);
+        valtemp.push_back(2.0 * x[idy]);
       }
     }
   }
@@ -116,29 +120,29 @@ void jac_f_ref(std::vector<T> x, std::vector<T> y, SparseMatrix& jac)
 }
 
 /// Function that computes the Jacobian via automatic differentiation
-template <typename T>
-__attribute__((noinline)) void jac_f(size_t N, T* input, SparseMatrix& jac)
+template <typename ScalarT>
+__attribute__((noinline)) void jac_f(size_t N, ScalarT* input, SparseMatrix& jac)
 {
-  std::vector<Triple<T>> triplets;
+  std::vector<Triple<ScalarT>> triplets;
   for (size_t i = 0; i < N; i++)
   {
-    T* output   = __enzyme_todense<T*>((void*) ident_load<T>, (void*) ident_store<T>, i);
-    T* d_output = __enzyme_todense<T*>((void*) sparse_load<T>, (void*) sparse_store<T>, i, &triplets);
+    ScalarT* output   = __enzyme_todense<ScalarT*>((void*) ident_load<ScalarT>, (void*) ident_store<ScalarT>, i);
+    ScalarT* d_output = __enzyme_todense<ScalarT*>((void*) sparse_load<ScalarT>, (void*) sparse_store<ScalarT>, i, &triplets);
 
-    __enzyme_fwddiff<void>((void*) f<T>,
+    __enzyme_fwddiff<void>((void*) f<ScalarT>,
                            enzyme_const,
                            N,
                            enzyme_dup,
                            input,
                            output,
                            enzyme_dupnoneed,
-                           (T*) 0x1,
+                           (ScalarT*) 0x1,
                            d_output);
   }
 
-  std::vector<size_t> ctemp{};
-  std::vector<size_t> rtemp{};
-  std::vector<T>      valtemp{};
+  std::vector<size_t>  ctemp{};
+  std::vector<size_t>  rtemp{};
+  std::vector<ScalarT> valtemp{};
   for (auto& tup : triplets)
   {
     rtemp.push_back(tup.row);
@@ -155,13 +159,13 @@ void check(SparseMatrix matrix_1, SparseMatrix matrix_2, int& fail)
   const auto [rcord_1, ccord_1, vals_1]                                                  = entries_1;
   std::tuple<std::vector<size_t>&, std::vector<size_t>&, std::vector<double>&> entries_2 = matrix_2.getEntries();
   const auto [rcord_2, ccord_2, vals_2]                                                  = entries_2;
-  for (int ind = 0; ind < vals_1.size(); ++ind)
+  for (size_t ind = 0; ind < vals_1.size(); ++ind)
   {
     if (rcord_1[ind] != rcord_2[ind])
       fail++;
     if (ccord_1[ind] != ccord_2[ind])
       fail++;
-    if (std::abs(vals_1[ind] - vals_2[ind]) > std::numeric_limits<double>::epsilon())
+    if (!GridKit::Testing::isEqual(vals_1[ind], vals_2[ind]))
       fail++;
   }
 }
@@ -177,7 +181,7 @@ int main()
 
   /// Input initialization
   double val = 0.0;
-  for (int i = 0; i < N; ++i)
+  for (size_t i = 0; i < N; ++i)
   {
     x[i]  = val;
     val  += 1.0;
diff --git a/examples/Enzyme/Standalone/EnzymeVector.cpp b/examples/Enzyme/Standalone/EnzymeVector.cpp
index 94a10f9a..4883fc56 100644
--- a/examples/Enzyme/Standalone/EnzymeVector.cpp
+++ b/examples/Enzyme/Standalone/EnzymeVector.cpp
@@ -1,8 +1,8 @@
 #include <iostream>
-#include <limits>
 #include <vector>
 
 #include <LinearAlgebra/DenseMatrix/DenseMatrix.hpp>
+#include <Utilities/Testing.hpp>
 
 /**
  * @brief Standalone example that computes the Jacobian of a vector-valued function
@@ -30,21 +30,25 @@ inline double dsquare_ref_scalar(double x)
 // Vector-valued function to differentiate
 void square(std::vector<double> x, std::vector<double>& y)
 {
-  for (int idx = 0; idx < x.size(); ++idx)
+  for (size_t idx = 0; idx < x.size(); ++idx)
   {
-    y[idx] = square_scalar(x[idx]);
+    y[idx] = 0.0;
+    for (size_t idy = 0; idy <= idx; idy++)
+    {
+      y[idx] += square_scalar(x[idy]);
+    }
   }
 }
 
 // Reference Jacobian
 void dsquare_ref(std::vector<double> x, std::vector<double> y, DenseMatrix& dy)
 {
-  for (int idy = 0; idy < y.size(); ++idy)
+  for (size_t idy = 0; idy < y.size(); ++idy)
   {
-    for (int idx = 0; idx < x.size(); ++idx)
+    for (size_t idx = 0; idx < x.size(); ++idx)
     {
-      if (idx == idy)
-        dy.setValue(idx, idy, dsquare_ref_scalar(x[idx]));
+      if (idy <= idx)
+        dy.setValue(idx, idy, dsquare_ref_scalar(x[idy]));
     }
   }
 }
@@ -54,10 +58,10 @@ void dsquare(std::vector<double> x, std::vector<double> y, DenseMatrix& dy)
 {
   std::vector<double> v(x.size());
   std::vector<double> d_y(y.size());
-  for (int idy = 0; idy < y.size(); ++idy)
+  for (size_t idy = 0; idy < y.size(); ++idy)
   {
     // Elementary vector for Jacobian-vector product
-    for (int idx = 0; idx < x.size(); ++idx)
+    for (size_t idx = 0; idx < x.size(); ++idx)
     {
       v[idx] = 0.0;
     }
@@ -67,7 +71,7 @@ void dsquare(std::vector<double> x, std::vector<double> y, DenseMatrix& dy)
     __enzyme_fwddiff((void*) square, enzyme_dup, x, v, enzyme_dupnoneed, y, &d_y);
 
     // Store result
-    for (int idx = 0; idx < x.size(); ++idx)
+    for (size_t idx = 0; idx < x.size(); ++idx)
     {
       dy.setValue(idx, idy, d_y[idx]);
     }
@@ -77,15 +81,15 @@ void dsquare(std::vector<double> x, std::vector<double> y, DenseMatrix& dy)
 int main()
 {
   // Vector and matrix declarations
-  constexpr int       N = 10;
+  constexpr size_t    N = 10;
   std::vector<double> x(N);
   std::vector<double> sq(N);
   DenseMatrix         dsq     = DenseMatrix(N, N);
   DenseMatrix         dsq_ref = DenseMatrix(N, N);
 
   // Random input values
-  srand(time(NULL));
-  for (int idx = 0; idx < x.size(); ++idx)
+  srand(static_cast<unsigned int>(time(NULL)));
+  for (size_t idx = 0; idx < x.size(); ++idx)
   {
     x[idx] = rand();
   }
@@ -102,11 +106,11 @@ int main()
   // Check
   int  fail    = 0;
   bool verbose = true;
-  for (int idy = 0; idy < sq.size(); ++idy)
+  for (size_t idy = 0; idy < sq.size(); ++idy)
   {
-    for (int idx = 0; idx < x.size(); ++idx)
+    for (size_t idx = 0; idx < x.size(); ++idx)
     {
-      if (std::abs(dsq.getValue(idx, idy) - dsq_ref.getValue(idx, idy)) > std::numeric_limits<double>::epsilon())
+      if (!GridKit::Testing::isEqual(dsq.getValue(idx, idy), dsq_ref.getValue(idx, idy)))
       {
         fail++;
         if (verbose)
diff --git a/src/AutomaticDifferentiation/Enzyme/SparseWrapper.hpp b/src/AutomaticDifferentiation/Enzyme/SparseWrapper.hpp
new file mode 100644
index 00000000..61ec7ee9
--- /dev/null
+++ b/src/AutomaticDifferentiation/Enzyme/SparseWrapper.hpp
@@ -0,0 +1,176 @@
+#pragma once
+
+#include <vector>
+
+/**
+ * @brief Enzyme constants for activity analysis
+ *
+ */
+extern int enzyme_dup;
+extern int enzyme_const;
+extern int enzyme_dupnoneed;
+
+namespace GridKit
+{
+  namespace Enzyme
+  {
+    /**
+     * @brief Residual wrapper around residual methods inside model classes
+     *
+     * @tparam ModelT - model type
+     * @tparam ScalarT - scalar data type
+     */
+    template <typename ModelT, typename ScalarT>
+    void residual_wrapper(ModelT* obj, ScalarT* y, ScalarT* f)
+    {
+      obj->evaluateResidualLocally(y, f);
+    }
+
+    /**
+     * @brief Enzyme fwddiff template
+     *
+     * @tparam T - return type
+     * @tparam ModelT - model type
+     */
+    template <typename T, typename... ModelT>
+    extern T __enzyme_fwddiff(void*, ModelT...) noexcept;
+
+    /**
+     * @brief Enzyme todense template
+     *
+     * @tparam T - return type
+     * @tparam ModelT - model type
+     */
+    template <typename T, typename... ModelT>
+    extern T __enzyme_todense(ModelT...) noexcept;
+
+    /**
+     * @brief Enzyme sparse storage in triplet format
+     *
+     * @tparam ScalarT - scalar data type
+     */
+    template <typename ScalarT>
+    struct Triple
+    {
+      size_t  row;
+      size_t  col;
+      ScalarT val;
+      Triple(Triple&&) = default;
+
+      Triple(size_t row, size_t col, ScalarT val)
+        : row(row),
+          col(col),
+          val(val)
+      {
+      }
+    };
+
+    /**
+     * @brief Enzyme sparse accumulation for float
+     *
+     */
+    [[maybe_unused]] __attribute__((enzyme_sparse_accumulate)) static void inner_storeflt(size_t row, size_t col, float val, std::vector<Triple<float>>& triplets)
+    {
+      triplets.emplace_back(row, col, val);
+    }
+
+    /**
+     * @brief Enzyme sparse accumulation for double
+     *
+     */
+    [[maybe_unused]] __attribute__((enzyme_sparse_accumulate)) static void inner_storedbl(size_t row, size_t col, double val, std::vector<Triple<double>>& triplets)
+    {
+      triplets.emplace_back(row, col, val);
+    }
+
+    /**
+     * @brief Enzyme sparse store
+     *
+     * @tparam ScalarT - scalar data type
+     */
+    template <typename ScalarT>
+    __attribute__((always_inline)) static void sparse_store(ScalarT val, size_t idx, size_t i, std::vector<Triple<ScalarT>>& triplets)
+    {
+      if (val == 0.0)
+        return;
+      idx /= sizeof(ScalarT);
+      if constexpr (sizeof(ScalarT) == 4)
+        inner_storeflt(idx, i, val, triplets);
+      else
+        inner_storedbl(idx, i, val, triplets);
+    }
+
+    /**
+     * @brief Enzyme sparse load
+     *
+     * @tparam ScalarT - scalar data type
+     */
+    template <typename ScalarT>
+    __attribute__((always_inline)) static ScalarT sparse_load(size_t, size_t, std::vector<Triple<ScalarT>>&)
+    {
+      return 0.0;
+    }
+
+    /**
+     * @brief Enzyme identity store
+     *
+     * @tparam ScalarT - scalar data type
+     */
+    template <typename ScalarT>
+    __attribute__((always_inline)) static void ident_store(ScalarT, size_t, size_t)
+    {
+      assert(0 && "should never load");
+    }
+
+    /**
+     * @brief Enzyme identity load
+     *
+     * @tparam ScalarT - scalar data type
+     */
+    template <typename ScalarT>
+    __attribute__((always_inline)) static ScalarT ident_load(size_t idx, size_t i)
+    {
+      idx /= sizeof(ScalarT);
+      return (ScalarT) (idx == i);
+    }
+
+    /**
+     * @brief Function that computes the Jacobian via automatic differentiation
+     *
+     * @tparam ModelT - model type
+     * @tparam ScalarT - scalar data type
+     * @tparam IdxT    - matrix index data type
+     */
+    template <typename ModelT, class ScalarT, typename IdxT>
+    __attribute__((noinline)) void EnzymeSparseModelJacobian(ModelT* model, size_t n, ScalarT* input, GridKit::LinearAlgebra::COO_Matrix<ScalarT, IdxT>& jac)
+    {
+      std::vector<Triple<ScalarT>> triplets;
+      for (size_t i = 0; i < n; i++)
+      {
+        ScalarT* output   = __enzyme_todense<ScalarT*>((void*) ident_load<ScalarT>, (void*) ident_store<ScalarT>, i);
+        ScalarT* d_output = __enzyme_todense<ScalarT*>((void*) sparse_load<ScalarT>, (void*) sparse_store<ScalarT>, i, &triplets);
+
+        __enzyme_fwddiff<void>((void*) residual_wrapper<ModelT, ScalarT>,
+                               enzyme_const,
+                               model,
+                               enzyme_dup,
+                               input,
+                               output,
+                               enzyme_dupnoneed,
+                               (ScalarT*) 0x1,
+                               d_output);
+      }
+
+      std::vector<IdxT>    ctemp{};
+      std::vector<IdxT>    rtemp{};
+      std::vector<ScalarT> valtemp{};
+      for (auto& tup : triplets)
+      {
+        rtemp.push_back(static_cast<IdxT>(tup.row));
+        ctemp.push_back(static_cast<IdxT>(tup.col));
+        valtemp.push_back(tup.val);
+      }
+      jac.setValues(rtemp, ctemp, valtemp);
+    }
+  } // namespace Enzyme
+} // namespace GridKit
diff --git a/src/Model/PhasorDynamics/Load/CMakeLists.txt b/src/Model/PhasorDynamics/Load/CMakeLists.txt
index 0f1659a7..909d04dc 100644
--- a/src/Model/PhasorDynamics/Load/CMakeLists.txt
+++ b/src/Model/PhasorDynamics/Load/CMakeLists.txt
@@ -1,12 +1,29 @@
- 
+
 
 # [[
 #  Author(s):
 #    - Cameron Rutherford <cameron.rutherford@pnnl.gov>
 # ]]
 
-gridkit_add_library(phasor_dynamics_load
+
+if(GRIDKIT_ENABLE_ENZYME)
+  gridkit_add_library(phasor_dynamics_load
+    SOURCES
+      LoadEnzyme.cpp
+    LINK_LIBRARIES
+      ClangEnzymeFlags
+    OUTPUT_NAME
+      gridkit_phasor_dynamics_load)
+else()
+  gridkit_add_library(phasor_dynamics_load
+    SOURCES
+      Load.cpp
+    OUTPUT_NAME
+      gridkit_phasor_dynamics_load)
+endif()
+
+gridkit_add_library(phasor_dynamics_load_dependency_tracking
   SOURCES
-    Load.cpp
+    LoadDependencyTracking.cpp
   OUTPUT_NAME
-    gridkit_phasor_dynamics_load)
+    gridkit_phasor_dynamics_load_dependency_tracking)
diff --git a/src/Model/PhasorDynamics/Load/Load.cpp b/src/Model/PhasorDynamics/Load/Load.cpp
index d3f9e016..bbc18790 100644
--- a/src/Model/PhasorDynamics/Load/Load.cpp
+++ b/src/Model/PhasorDynamics/Load/Load.cpp
@@ -1,187 +1,29 @@
-
-#include "Load.hpp"
-
-#include <cmath>
-#include <iostream>
-
-#include <Model/PhasorDynamics/Bus/Bus.hpp>
-#include <Model/PhasorDynamics/Load/LoadData.hpp>
-
-namespace GridKit
-{
-  namespace PhasorDynamics
-  {
-    /*!
-     * @brief Constructor for a pi-model load
-     *
-     * Arguments passed to ModelEvaluatorImpl:
-     * - Number of equations = 0
-     * - Number of independent variables = 0
-     * - Number of quadratures = 0
-     * - Number of optimization parameters = 0
-     */
-
-    template <class ScalarT, typename IdxT>
-    Load<ScalarT, IdxT>::Load(bus_type* bus)
-      : bus_(bus)
-    {
-      size_ = 0;
-    }
-
-    template <class ScalarT, typename IdxT>
-    Load<ScalarT, IdxT>::Load(bus_type* bus,
-                              real_type R,
-                              real_type X)
-      : bus_(bus),
-        R_(R),
-        X_(X)
-    {
-    }
-
-    template <class ScalarT, typename IdxT>
-    Load<ScalarT, IdxT>::Load(bus_type*        bus,
-                              model_data_type& data)
-      : bus_(bus),
-        R_(data.R),
-        X_(data.X)
-    {
-    }
-
-    template <class ScalarT, typename IdxT>
-    Load<ScalarT, IdxT>::Load(bus_type* bus, IdxT component_id)
-      : bus_(bus)
-    {
-      size_         = 0;
-      component_id_ = component_id;
-    }
-
-    template <class ScalarT, typename IdxT>
-    Load<ScalarT, IdxT>::~Load()
-    {
-      // std::cout << "Destroy Load..." << std::endl;
-    }
-
-    /*!
-     * @brief allocate method computes sparsity pattern of the Jacobian.
-     */
-    template <class ScalarT, typename IdxT>
-    int Load<ScalarT, IdxT>::allocate()
-    {
-      // std::cout << "Allocate Load..." << std::endl;
-      return 0;
-    }
-
-    /**
-     * Initialization of the load model
-     *
-     */
-    template <class ScalarT, typename IdxT>
-    int Load<ScalarT, IdxT>::initialize()
-    {
-      return 0;
-    }
-
-    /**
-     * \brief Identify differential variables.
-     */
-    template <class ScalarT, typename IdxT>
-    int Load<ScalarT, IdxT>::tagDifferentiable()
-    {
-      return 0;
-    }
-
-    /**
-     * \brief Residual contribution of the load is pushed to the bus.
-     *
-     */
-    template <class ScalarT, typename IdxT>
-    int Load<ScalarT, IdxT>::evaluateResidual()
-    {
-      real_type b = -X_ / (R_ * R_ + X_ * X_);
-      real_type g = R_ / (R_ * R_ + X_ * X_);
-
-      Ir() += -g * Vr() + b * Vi();
-      Ii() += -b * Vr() - g * Vi();
-
-      return 0;
-    }
-
-    /**
-     * @brief Jacobian evaluation not implemented yet
-     *
-     * @tparam ScalarT - scalar data type
-     * @tparam IdxT    - matrix index data type
-     * @return int - error code, 0 = success
-     */
-    template <class ScalarT, typename IdxT>
-    int Load<ScalarT, IdxT>::evaluateJacobian()
-    {
-      std::cout << "Evaluate Jacobian for Load..." << std::endl;
-      std::cout << "Jacobian evaluation not implemented!" << std::endl;
-      return 0;
-    }
-
-    /**
-     * @brief Integrand (objective) evaluation not implemented yet
-     *
-     * @tparam ScalarT - scalar data type
-     * @tparam IdxT    - matrix index data type
-     * @return int - error code, 0 = success
-     */
-    template <class ScalarT, typename IdxT>
-    int Load<ScalarT, IdxT>::evaluateIntegrand()
-    {
-      // std::cout << "Evaluate Integrand for Load..." << std::endl;
-      return 0;
-    }
-
-    /**
-     * @brief Adjoint initialization not implemented yet
-     *
-     * @tparam ScalarT - scalar data type
-     * @tparam IdxT    - matrix index data type
-     * @return int - error code, 0 = success
-     */
-    template <class ScalarT, typename IdxT>
-    int Load<ScalarT, IdxT>::initializeAdjoint()
-    {
-      // std::cout << "Initialize adjoint for Load..." << std::endl;
-      return 0;
-    }
-
-    /**
-     * @brief Adjoint residual evaluation not implemented yet
-     *
-     * @tparam ScalarT - scalar data type
-     * @tparam IdxT    - matrix index data type
-     * @return int - error code, 0 = success
-     */
-    template <class ScalarT, typename IdxT>
-    int Load<ScalarT, IdxT>::evaluateAdjointResidual()
-    {
-      // std::cout << "Evaluate adjoint residual for Load..." << std::endl;
-      return 0;
-    }
-
-    /**
-     * @brief Adjoint integrand (objective) evaluation not implemented yet
-     *
-     * @tparam ScalarT - scalar data type
-     * @tparam IdxT    - matrix index data type
-     * @return int - error code, 0 = success
-     */
-    template <class ScalarT, typename IdxT>
-    int Load<ScalarT, IdxT>::evaluateAdjointIntegrand()
-    {
-      // std::cout << "Evaluate adjoint Integrand for Load..." << std::endl;
-      return 0;
-    }
-
-    // Available template instantiations
-    template class Load<double, long int>;
-    template class Load<double, size_t>;
-    template class Load<DependencyTracking::Variable, long int>;
-    template class Load<DependencyTracking::Variable, size_t>;
-
-  } // namespace PhasorDynamics
-} // namespace GridKit
+
+#include "LoadImpl.hpp"
+
+namespace GridKit
+{
+  namespace PhasorDynamics
+  {
+    /**
+     * @brief Jacobian evaluation not implemented
+     *
+     * @tparam ScalarT - scalar data type
+     * @tparam IdxT    - matrix index data type
+     * @return int - error code, 0 = success
+     */
+    template <class ScalarT, typename IdxT>
+    int Load<ScalarT, IdxT>::evaluateJacobian()
+    {
+      std::cout << "Evaluate Jacobian for Load..." << std::endl;
+      std::cout << "Jacobian evaluation is not implemented!" << std::endl;
+
+      return 0;
+    }
+
+    // Available template instantiations
+    template class Load<double, long int>;
+    template class Load<double, size_t>;
+
+  } // namespace PhasorDynamics
+} // namespace GridKit
diff --git a/src/Model/PhasorDynamics/Load/Load.hpp b/src/Model/PhasorDynamics/Load/Load.hpp
index 6fc996c8..ee20f428 100644
--- a/src/Model/PhasorDynamics/Load/Load.hpp
+++ b/src/Model/PhasorDynamics/Load/Load.hpp
@@ -36,6 +36,7 @@ namespace GridKit
       using Component<ScalarT, IdxT>::yp_;
       using Component<ScalarT, IdxT>::tag_;
       using Component<ScalarT, IdxT>::f_;
+      using Component<ScalarT, IdxT>::J_;
       using Component<ScalarT, IdxT>::component_id_;
 
       using real_type       = typename Component<ScalarT, IdxT>::real_type;
@@ -98,6 +99,9 @@ namespace GridKit
         return bus_->Ii();
       }
 
+    public:
+      int evaluateResidualLocally(ScalarT*, ScalarT*);
+
     private:
       bus_type* bus_{nullptr};
       real_type R_{0.1};
diff --git a/src/Model/PhasorDynamics/Load/LoadDependencyTracking.cpp b/src/Model/PhasorDynamics/Load/LoadDependencyTracking.cpp
new file mode 100644
index 00000000..8219abca
--- /dev/null
+++ b/src/Model/PhasorDynamics/Load/LoadDependencyTracking.cpp
@@ -0,0 +1,29 @@
+
+#include "LoadImpl.hpp"
+
+namespace GridKit
+{
+  namespace PhasorDynamics
+  {
+    /**
+     * @brief Jacobian evaluation not implemented
+     *
+     * @tparam ScalarT - scalar data type
+     * @tparam IdxT    - matrix index data type
+     * @return int - error code, 0 = success
+     */
+    template <class ScalarT, typename IdxT>
+    int Load<ScalarT, IdxT>::evaluateJacobian()
+    {
+      std::cout << "Evaluate Jacobian for Load..." << std::endl;
+      std::cout << "Jacobian evaluation is not implemented!" << std::endl;
+
+      return 0;
+    }
+
+    // Available template instantiations
+    template class Load<DependencyTracking::Variable, long int>;
+    template class Load<DependencyTracking::Variable, size_t>;
+
+  } // namespace PhasorDynamics
+} // namespace GridKit
diff --git a/src/Model/PhasorDynamics/Load/LoadEnzyme.cpp b/src/Model/PhasorDynamics/Load/LoadEnzyme.cpp
new file mode 100644
index 00000000..1410290f
--- /dev/null
+++ b/src/Model/PhasorDynamics/Load/LoadEnzyme.cpp
@@ -0,0 +1,39 @@
+
+#include "LoadImpl.hpp"
+#include <AutomaticDifferentiation/Enzyme/SparseWrapper.hpp>
+
+namespace GridKit
+{
+  namespace PhasorDynamics
+  {
+    /**
+     * @brief Jacobian evaluation experimental
+     *
+     * @tparam ScalarT - scalar data type
+     * @tparam IdxT    - matrix index data type
+     * @return int - error code, 0 = success
+     */
+    template <class ScalarT, typename IdxT>
+    int Load<ScalarT, IdxT>::evaluateJacobian()
+    {
+      std::cout << "Evaluate Jacobian for Load..." << std::endl;
+      std::cout << "Jacobian evaluation is experimental!" << std::endl;
+
+      std::vector<ScalarT> y(2);
+      std::vector<ScalarT> f(2);
+      y[0] = Vr();
+      y[1] = Vi();
+      /// Setting J_ via Enzyme works, though J_ had not been initialized within the model.
+      /// This is because the COO_Matrix class is very permissive.
+      /// Having the currents as model variables and allocating J_ accordingly will be helpful.
+      GridKit::Enzyme::EnzymeSparseModelJacobian<Load<ScalarT, IdxT>, ScalarT, IdxT>(this, f.size(), y.data(), J_);
+
+      return 0;
+    }
+
+    // Available template instantiations
+    template class Load<double, long int>;
+    template class Load<double, size_t>;
+
+  } // namespace PhasorDynamics
+} // namespace GridKit
diff --git a/src/Model/PhasorDynamics/Load/LoadImpl.hpp b/src/Model/PhasorDynamics/Load/LoadImpl.hpp
new file mode 100644
index 00000000..e8cfbd2b
--- /dev/null
+++ b/src/Model/PhasorDynamics/Load/LoadImpl.hpp
@@ -0,0 +1,183 @@
+
+#include <cmath>
+#include <iostream>
+
+#include "Load.hpp"
+#include <Model/PhasorDynamics/Bus/Bus.hpp>
+#include <Model/PhasorDynamics/Load/LoadData.hpp>
+
+namespace GridKit
+{
+  namespace PhasorDynamics
+  {
+    /*!
+     * @brief Constructor for a pi-model load
+     *
+     * Arguments passed to ModelEvaluatorImpl:
+     * - Number of equations = 0
+     * - Number of independent variables = 0
+     * - Number of quadratures = 0
+     * - Number of optimization parameters = 0
+     */
+
+    template <class ScalarT, typename IdxT>
+    Load<ScalarT, IdxT>::Load(bus_type* bus)
+      : bus_(bus)
+    {
+      size_ = 0;
+    }
+
+    template <class ScalarT, typename IdxT>
+    Load<ScalarT, IdxT>::Load(bus_type* bus,
+                              real_type R,
+                              real_type X)
+      : bus_(bus),
+        R_(R),
+        X_(X)
+    {
+    }
+
+    template <class ScalarT, typename IdxT>
+    Load<ScalarT, IdxT>::Load(bus_type*        bus,
+                              model_data_type& data)
+      : bus_(bus),
+        R_(data.R),
+        X_(data.X)
+    {
+    }
+
+    template <class ScalarT, typename IdxT>
+    Load<ScalarT, IdxT>::Load(bus_type* bus, IdxT component_id)
+      : bus_(bus)
+    {
+      size_         = 0;
+      component_id_ = component_id;
+    }
+
+    template <class ScalarT, typename IdxT>
+    Load<ScalarT, IdxT>::~Load()
+    {
+      // std::cout << "Destroy Load..." << std::endl;
+    }
+
+    /*!
+     * @brief allocate method computes sparsity pattern of the Jacobian.
+     */
+    template <class ScalarT, typename IdxT>
+    int Load<ScalarT, IdxT>::allocate()
+    {
+      // std::cout << "Allocate Load..." << std::endl;
+      return 0;
+    }
+
+    /**
+     * Initialization of the load model
+     *
+     */
+    template <class ScalarT, typename IdxT>
+    int Load<ScalarT, IdxT>::initialize()
+    {
+      return 0;
+    }
+
+    /**
+     * \brief Identify differential variables.
+     */
+    template <class ScalarT, typename IdxT>
+    int Load<ScalarT, IdxT>::tagDifferentiable()
+    {
+      return 0;
+    }
+
+    /**
+     * @brief Residual contribution computed locally
+     *
+     */
+    template <class ScalarT, typename IdxT>
+    int Load<ScalarT, IdxT>::evaluateResidualLocally(ScalarT* y, ScalarT* f)
+    {
+      real_type b = -X_ / (R_ * R_ + X_ * X_);
+      real_type g = R_ / (R_ * R_ + X_ * X_);
+
+      f[0] = -g * y[0] + b * y[1];
+      f[1] = -b * y[0] - g * y[1];
+
+      return 0;
+    }
+
+    /**
+     * @brief Residual contribution of the load is pushed to the bus.
+     *
+     */
+    template <class ScalarT, typename IdxT>
+    int Load<ScalarT, IdxT>::evaluateResidual()
+    {
+      std::vector<ScalarT> y(2);
+      std::vector<ScalarT> f(2);
+      y[0] = Vr();
+      y[1] = Vi();
+      evaluateResidualLocally(y.data(), f.data());
+      Ir() += f[0];
+      Ii() += f[1];
+
+      return 0;
+    }
+
+    /**
+     * @brief Integrand (objective) evaluation not implemented yet
+     *
+     * @tparam ScalarT - scalar data type
+     * @tparam IdxT    - matrix index data type
+     * @return int - error code, 0 = success
+     */
+    template <class ScalarT, typename IdxT>
+    int Load<ScalarT, IdxT>::evaluateIntegrand()
+    {
+      // std::cout << "Evaluate Integrand for Load..." << std::endl;
+      return 0;
+    }
+
+    /**
+     * @brief Adjoint initialization not implemented yet
+     *
+     * @tparam ScalarT - scalar data type
+     * @tparam IdxT    - matrix index data type
+     * @return int - error code, 0 = success
+     */
+    template <class ScalarT, typename IdxT>
+    int Load<ScalarT, IdxT>::initializeAdjoint()
+    {
+      // std::cout << "Initialize adjoint for Load..." << std::endl;
+      return 0;
+    }
+
+    /**
+     * @brief Adjoint residual evaluation not implemented yet
+     *
+     * @tparam ScalarT - scalar data type
+     * @tparam IdxT    - matrix index data type
+     * @return int - error code, 0 = success
+     */
+    template <class ScalarT, typename IdxT>
+    int Load<ScalarT, IdxT>::evaluateAdjointResidual()
+    {
+      // std::cout << "Evaluate adjoint residual for Load..." << std::endl;
+      return 0;
+    }
+
+    /**
+     * @brief Adjoint integrand (objective) evaluation not implemented yet
+     *
+     * @tparam ScalarT - scalar data type
+     * @tparam IdxT    - matrix index data type
+     * @return int - error code, 0 = success
+     */
+    template <class ScalarT, typename IdxT>
+    int Load<ScalarT, IdxT>::evaluateAdjointIntegrand()
+    {
+      // std::cout << "Evaluate adjoint Integrand for Load..." << std::endl;
+      return 0;
+    }
+
+  } // namespace PhasorDynamics
+} // namespace GridKit
diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt
index dfe3c93d..c8b0a536 100644
--- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt
+++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt
@@ -13,6 +13,7 @@ target_link_libraries(test_phasor_branch GRIDKIT::phasor_dynamics_branch
                                 
 add_executable(test_phasor_load runLoadTests.cpp)
 target_link_libraries(test_phasor_load GRIDKIT::phasor_dynamics_load
+                                       GRIDKIT::phasor_dynamics_load_dependency_tracking
                                        GRIDKIT::phasor_dynamics_bus)
 
 add_executable(test_phasor_genrou runGenrouTests.cpp)
diff --git a/tests/UnitTests/PhasorDynamics/LoadTests.hpp b/tests/UnitTests/PhasorDynamics/LoadTests.hpp
index f9104336..081db141 100644
--- a/tests/UnitTests/PhasorDynamics/LoadTests.hpp
+++ b/tests/UnitTests/PhasorDynamics/LoadTests.hpp
@@ -100,6 +100,36 @@ namespace GridKit
         return success.report(__func__);
       }
 
+#ifdef GRIDKIT_ENABLE_ENZYME
+      TestOutcome enzyme_jacobian()
+      {
+        TestStatus success = true;
+
+        real_type R{2.0}; ///< Load resistance
+        real_type X{4.0}; ///< Load reactance
+
+        ScalarT Vr{10.0}; ///< Bus real voltage
+        ScalarT Vi{20.0}; ///< Bus imaginary voltage
+
+        PhasorDynamics::BusInfinite<ScalarT, IdxT> bus(Vr, Vi);
+
+        PhasorDynamics::Load<ScalarT, IdxT> load(&bus, R, X);
+        load.evaluateJacobian();
+        GridKit::LinearAlgebra::COO_Matrix<ScalarT, IdxT> model_jacobian = load.getJacobian();
+        model_jacobian.printMatrix("Model Jacobian");
+
+        /// Compare model Jacobian wih dependencies computed analytically
+        std::vector<DependencyTracking::Variable::DependencyMap> ref                = analyticalJacobian(R, X);
+        std::vector<DependencyTracking::Variable::DependencyMap> model_dependencies = mapFromCOO(model_jacobian);
+        for (size_t i = 0; i < ref.size(); ++i)
+        {
+          success *= (GridKit::Testing::isEqual(model_dependencies[i], ref[i]));
+        }
+
+        return success.report(__func__);
+      }
+#endif
+
     private:
       std::vector<DependencyTracking::Variable::DependencyMap> analyticalJacobian(const real_type R,
                                                                                   const real_type X)
@@ -119,6 +149,24 @@ namespace GridKit
 
         return dependencies;
       }
+
+      std::vector<DependencyTracking::Variable::DependencyMap> mapFromCOO(GridKit::LinearAlgebra::COO_Matrix<ScalarT, IdxT> matrix)
+      {
+        std::tuple<std::vector<IdxT>&, std::vector<IdxT>&, std::vector<ScalarT>&> matrix_entries = matrix.getEntries();
+        const auto [rows, columns, values]                                                       = matrix_entries;
+
+        std::tuple<IdxT, IdxT> matrix_dimensions = matrix.getDimensions();
+        const auto [n_rows, n_columns]           = matrix_dimensions;
+
+        std::vector<DependencyTracking::Variable::DependencyMap> dependencies(n_rows);
+
+        for (IdxT i = 0; i < rows.size(); ++i)
+        {
+          dependencies[rows[i]].insert(std::make_pair(columns[i], values[i]));
+        }
+
+        return dependencies;
+      }
     };
 
   } // namespace Testing
diff --git a/tests/UnitTests/PhasorDynamics/runLoadTests.cpp b/tests/UnitTests/PhasorDynamics/runLoadTests.cpp
index 41426a82..1dac30e2 100644
--- a/tests/UnitTests/PhasorDynamics/runLoadTests.cpp
+++ b/tests/UnitTests/PhasorDynamics/runLoadTests.cpp
@@ -11,6 +11,9 @@ int main()
   result += test.constructor();
   result += test.residual();
   result += test.jacobian();
+#ifdef GRIDKIT_ENABLE_ENZYME
+  result += test.enzyme_jacobian();
+#endif
 
   return result.summary();
 }