diff --git a/cmake/FindEnzyme.cmake b/cmake/FindEnzyme.cmake
index bc1ebbdab..14420797e 100644
--- a/cmake/FindEnzyme.cmake
+++ b/cmake/FindEnzyme.cmake
@@ -101,7 +101,7 @@ message(STATUS "opt: ${GRIDKIT_OPT}")
 macro(enzyme_add_executable)
   set(options)
   set(oneValueArgs NAME)
-  set(multiValueArgs SOURCES LINK_LIBRARIES)
+  set(multiValueArgs SOURCES LINK_LIBRARIES INCLUDE_DIRECTORIES)
   cmake_parse_arguments(enzyme_add_executable "${options}" "${oneValueArgs}"
     "${multiValueArgs}" ${ARGN})
 
@@ -111,14 +111,32 @@ macro(enzyme_add_executable)
   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}
+      DEPENDS ${PHASE0} 
       OUTPUT ${PHASE1}
-      COMMAND ${CMAKE_CXX_COMPILER} -flto -c ${PHASE0} -O2 -fno-vectorize -ffast-math -fno-unroll-loops -o ${PHASE1}
+      COMMAND ${CMAKE_CXX_COMPILER} -flto -c ${PHASE0} ${INCLUDE_COMPILER_LIST} -O2 -fno-vectorize -ffast-math -fno-unroll-loops -o ${PHASE1}
       COMMENT "Compiling ${SRC} to object file for target ${enzyme_add_executable_NAME}"
       )
     set(OBJS "${OBJS} ${PHASE1}")
@@ -148,9 +166,9 @@ macro(enzyme_add_executable)
     )
 
   add_custom_command(
-    DEPENDS ${PHASE4}
+    DEPENDS ${PHASE4} ${enzyme_add_executable_LINK_LIBRARIES}
     OUTPUT ${PHASE5}
-    COMMAND ${CMAKE_CXX_COMPILER} ${PHASE4} -o ${PHASE5}
+    COMMAND ${CMAKE_CXX_COMPILER} ${LINKER_FLAGS} ${PHASE4} -o ${PHASE5}
     )
 
   add_custom_target(
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 5b75c689b..100de2338 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -56,7 +56,7 @@
 #
 
 add_subdirectory(MatPowerTesting)
-add_subdirectory(SparseTest)
+add_subdirectory(LinearAlgebra)
 add_subdirectory(DistributedGeneratorTest)
 
 if(TARGET SUNDIALS::kinsol)
diff --git a/examples/DistributedGeneratorTest/CMakeLists.txt b/examples/DistributedGeneratorTest/CMakeLists.txt
index 41b124cb8..6ed544e93 100644
--- a/examples/DistributedGeneratorTest/CMakeLists.txt
+++ b/examples/DistributedGeneratorTest/CMakeLists.txt
@@ -3,9 +3,9 @@
 
 
 add_executable(dgtest DGTest.cpp)
-target_link_libraries(dgtest GRIDKIT::powerelec_disgen
-                                GRIDKIT::powerelec_mircoline
-                                GRIDKIT::powerelec_microload
+target_link_libraries(dgtest GRIDKIT::power_elec_disgen
+                                GRIDKIT::power_elec_microline
+                                GRIDKIT::power_elec_microload
                                 GRIDKIT::solvers_dyn)
                                 
 add_test(NAME DistributedGeneratorTest    COMMAND $<TARGET_FILE:dgtest>)
diff --git a/examples/Enzyme/CMakeLists.txt b/examples/Enzyme/CMakeLists.txt
index 794ec850c..f7ed45025 100644
--- a/examples/Enzyme/CMakeLists.txt
+++ b/examples/Enzyme/CMakeLists.txt
@@ -1,2 +1,3 @@
 add_subdirectory(Standalone)
 add_subdirectory(Library)
+add_subdirectory(PowerElectronics)
diff --git a/examples/Enzyme/Library/CMakeLists.txt b/examples/Enzyme/Library/CMakeLists.txt
index 78566c983..6dec74892 100644
--- a/examples/Enzyme/Library/CMakeLists.txt
+++ b/examples/Enzyme/Library/CMakeLists.txt
@@ -1,4 +1,2 @@
-enzyme_add_executable(NAME EnzymeLibCheck SOURCES main.cpp library.cpp)
-#install(TARGETS ${CMAKE_CURRENT_BINARY_DIR}/EnzymeLibCheck DESTINATION bin)
-
-add_test(NAME "EnzymeLibCheck" COMMAND ${CMAKE_CURRENT_BINARY_DIR}/EnzymeLibCheck)
+add_subdirectory(Scalar)
+add_subdirectory(Vector)
diff --git a/examples/Enzyme/Library/Scalar/CMakeLists.txt b/examples/Enzyme/Library/Scalar/CMakeLists.txt
new file mode 100644
index 000000000..6dff430a6
--- /dev/null
+++ b/examples/Enzyme/Library/Scalar/CMakeLists.txt
@@ -0,0 +1,6 @@
+enzyme_add_executable(
+  NAME EnzymeLibScalarCheck 
+  SOURCES EnzymeScalar.cpp ScalarModel.cpp
+)
+
+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
new file mode 100644
index 000000000..abc1a9d07
--- /dev/null
+++ b/examples/Enzyme/Library/Scalar/EnzymeScalar.cpp
@@ -0,0 +1,32 @@
+#include <iostream>
+#include <limits>
+#include "ScalarModel.hpp"
+
+/**
+ * @brief Example that computes the derivative of a library function 
+ * (implemented as the member function of a class and operating directly on class members)
+ * by automatic differentiation via Enzyme.
+ *
+ * TODO: Convert this into a unit test.
+ */
+
+int main() 
+{
+    int fail = 0;
+    ScalarModel scalar_model;
+    double var = 5.0;
+    scalar_model.setVariable(var);
+    scalar_model.evalFunction();
+    scalar_model.evalDerivative();
+    double sq = scalar_model.getFunctionValue();
+    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())
+    {
+        fail++;
+        std::cout << "Result incorrect\n";
+    }
+    std::cout << "Status: " << fail << "\n";
+    return fail;
+}
diff --git a/examples/Enzyme/Library/Scalar/EnzymeWrapper.hpp b/examples/Enzyme/Library/Scalar/EnzymeWrapper.hpp
new file mode 100644
index 000000000..8e659dea4
--- /dev/null
+++ b/examples/Enzyme/Library/Scalar/EnzymeWrapper.hpp
@@ -0,0 +1,16 @@
+#pragma once
+
+int enzyme_dup;
+int enzyme_dupnoneed;
+int enzyme_out;
+int enzyme_const;
+
+template <typename return_type, typename ... T>
+return_type __enzyme_fwddiff(return_type*, int, T* ... );
+
+template <typename return_type, typename T>
+return_type wrapper(T* obj) 
+{
+    obj->evalFunction();
+    return obj->getFunctionValue();
+}
diff --git a/examples/Enzyme/Library/Scalar/ScalarModel.cpp b/examples/Enzyme/Library/Scalar/ScalarModel.cpp
new file mode 100644
index 000000000..06fb9c1b6
--- /dev/null
+++ b/examples/Enzyme/Library/Scalar/ScalarModel.cpp
@@ -0,0 +1,49 @@
+#include "EnzymeWrapper.hpp"
+#include "ScalarModel.hpp"
+#include <iostream>
+
+ScalarModel::ScalarModel() 
+{
+}
+
+inline
+double ScalarModel::square(double x) 
+{
+    return x * x;
+}
+
+void ScalarModel::setVariable(double x) 
+{
+    x_ = x;
+}
+
+void ScalarModel::evalFunction() 
+{
+    f_ = square(x_);
+}
+
+void ScalarModel::evalDerivative() 
+{
+    ScalarModel d_scalar_model;
+    d_scalar_model.setVariable(1.0);
+    df_dx_ = __enzyme_fwddiff<double, ScalarModel>((double*)wrapper<double, ScalarModel>, enzyme_dup, this, &d_scalar_model);
+}
+
+double ScalarModel::getVariable() const 
+{
+    return x_;
+}
+
+double ScalarModel::getFunctionValue() const 
+{
+    return f_;
+}
+
+double ScalarModel::getDerivativeValue() const 
+{
+    return df_dx_;
+}
+
+ScalarModel::~ScalarModel() 
+{
+}
diff --git a/examples/Enzyme/Library/Scalar/ScalarModel.hpp b/examples/Enzyme/Library/Scalar/ScalarModel.hpp
new file mode 100644
index 000000000..2632bf59d
--- /dev/null
+++ b/examples/Enzyme/Library/Scalar/ScalarModel.hpp
@@ -0,0 +1,22 @@
+#pragma once
+
+/**
+ * @brief Class providing methods to evaluate a function and its derivative.
+ * This is used to test automatic differentiation.
+ */
+class ScalarModel 
+{
+private:
+    double x_, f_, df_dx_;
+    inline double square(double);
+
+public:
+    ScalarModel();
+    void setVariable(double);
+    void evalFunction();
+    void evalDerivative();
+    double getVariable() const;
+    double getFunctionValue() const;
+    double getDerivativeValue() const;
+    ~ScalarModel(); 
+};
diff --git a/examples/Enzyme/Library/Vector/CMakeLists.txt b/examples/Enzyme/Library/Vector/CMakeLists.txt
new file mode 100644
index 000000000..eda952de8
--- /dev/null
+++ b/examples/Enzyme/Library/Vector/CMakeLists.txt
@@ -0,0 +1,7 @@
+enzyme_add_executable(
+  NAME EnzymeLibVectorCheck 
+  SOURCES EnzymeVector.cpp VectorModel.cpp
+  LINK_LIBRARIES GRIDKIT::DenseMatrix
+)
+
+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
new file mode 100644
index 000000000..f18ef3e2a
--- /dev/null
+++ b/examples/Enzyme/Library/Vector/EnzymeVector.cpp
@@ -0,0 +1,88 @@
+#include <iostream>
+#include <limits>
+#include "VectorModel.hpp"
+
+/**
+ * @brief Example that computes the Jacobian of a vector-valued residual 
+ * (implemented as the member function of a class and operating directly on class members)
+ * by automatic differentiation via Enzyme.
+ *
+ * TODO: Convert this into a unit test.
+ */
+
+inline
+double dsquare_ref_scalar(double x) 
+{
+    return 2.0 * x;
+}
+
+// Reference Jacobian
+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 (int idx = 0; idx < x.size(); ++idx)
+        {
+            if (idx == idy) 
+                jac.setValue(idx, idy, dsquare_ref_scalar(x[idx]));
+        }
+    }
+    return jac;
+}
+
+int main() 
+{
+    // Size and variable declarations
+    constexpr int n = 10;
+    std::vector<double> var(n);
+
+    // Random input values
+    srand(time(NULL));
+    for (int idx = 0; idx < var.size(); ++idx)
+    {
+        var[idx] = rand();
+    }
+
+    // Model
+    VectorModel* vector_model = new VectorModel(n);
+    vector_model->setVariable(var);
+    vector_model->evalResidual();
+    vector_model->evalJacobian();
+    std::vector<double> var_temp = vector_model->getVariable();
+    std::vector<double> res = vector_model->getResidual();
+    DenseMatrix jac = vector_model->getJacobian();
+
+    // Reference Jacobian
+    DenseMatrix jac_ref = dsquare_ref(var, res);
+  
+    // Check
+    int fail = 0;
+    bool verbose = true;
+    for (int idy = 0; idy < res.size(); ++idy)
+    {
+        for (int idx = 0; idx < var.size(); ++idx)
+        {
+            if (std::abs(jac.getValue(idx, idy) - jac_ref.getValue(idx, idy)) > std::numeric_limits<double>::epsilon())
+            {
+                fail++;
+                if (verbose)
+                {
+                    std::cout << "Result incorrect at line = " << idy << ", column = " << idx << "\n";
+                    std::cout << "x = " << var_temp[idx] << ", x^2 = " << res[idx] << ", d(x^2)/dx = " << jac.getValue(idx, idy) << "\n"; 
+                }
+            }
+        }
+    }
+    if (verbose)
+    {
+        jac.printMatrix("Autodiff Jacobian");
+        jac_ref.printMatrix("Reference Jacobian");
+    }
+    std::cout << "Status: " << fail << "\n";
+
+    // Cleanup
+    delete vector_model;
+
+    return fail;
+}
diff --git a/examples/Enzyme/Library/Vector/EnzymeWrapper.hpp b/examples/Enzyme/Library/Vector/EnzymeWrapper.hpp
new file mode 100644
index 000000000..3306ab957
--- /dev/null
+++ b/examples/Enzyme/Library/Vector/EnzymeWrapper.hpp
@@ -0,0 +1,18 @@
+#pragma once
+
+#include <vector>
+
+int enzyme_dup;
+int enzyme_dupnoneed;
+int enzyme_out;
+int enzyme_const;
+
+template <typename T>
+std::vector<double> __enzyme_fwddiff(std::vector<double>*, int, T*, T*);
+
+template <typename T>
+std::vector<double> wrapper(T* obj) 
+{
+    obj->evalResidual();
+    return obj->getResidual();
+}
diff --git a/examples/Enzyme/Library/Vector/VectorModel.cpp b/examples/Enzyme/Library/Vector/VectorModel.cpp
new file mode 100644
index 000000000..33f312ccd
--- /dev/null
+++ b/examples/Enzyme/Library/Vector/VectorModel.cpp
@@ -0,0 +1,84 @@
+#include "EnzymeWrapper.hpp"
+#include "VectorModel.hpp"
+#include <iostream>
+
+VectorModel::VectorModel(int n) :
+   x_(n),
+   f_(n),
+   df_dx_(n, n)
+{
+}
+
+inline
+double VectorModel::square_scalar(double x) 
+{
+    return x * x;
+}
+
+void VectorModel::square(std::vector<double>& x, std::vector<double>& y) 
+{
+    for (int idx = 0; idx < x.size(); ++idx)
+    {
+        y[idx] = this->square_scalar(x[idx]);
+    }
+}
+
+void VectorModel::setVariable(std::vector<double> x) 
+{
+    for (int idx = 0; idx < x.size(); ++idx)
+    {
+        x_[idx] = x[idx];
+    }
+}
+
+void VectorModel::evalResidual() 
+{
+    square(x_, f_);
+}
+
+void VectorModel::evalJacobian() 
+{
+    const int n = x_.size();
+    std::vector<double> v(n);
+    VectorModel d_vector_model(n);
+    for (int idy = 0; idy < n; ++idy)
+    {
+        // Elementary vector for Jacobian-vector product
+        for (int idx = 0; idx < n; ++idx)
+        {
+            v[idx] = 0.0;
+        }
+        v[idy] = 1.0;
+        d_vector_model.setVariable(v);
+  
+        // Autodiff
+        std::vector<double> d_res = __enzyme_fwddiff<VectorModel>(
+                                      (std::vector<double>*)wrapper<VectorModel>, 
+                                      enzyme_dup, this, &d_vector_model);
+  
+        // Store result
+        for (int idx = 0; idx < n; ++idx)
+        {
+            df_dx_.setValue(idx, idy, d_res[idx]);
+        }
+    }
+}
+
+std::vector<double>& VectorModel::getVariable() 
+{
+    return x_;
+}
+
+std::vector<double>& VectorModel::getResidual() 
+{
+    return f_;
+}
+
+DenseMatrix& VectorModel::getJacobian() 
+{
+    return df_dx_;
+}
+
+VectorModel::~VectorModel() 
+{
+}
diff --git a/examples/Enzyme/Library/Vector/VectorModel.hpp b/examples/Enzyme/Library/Vector/VectorModel.hpp
new file mode 100644
index 000000000..5ff5d87eb
--- /dev/null
+++ b/examples/Enzyme/Library/Vector/VectorModel.hpp
@@ -0,0 +1,28 @@
+#pragma once
+
+#include <LinearAlgebra/DenseMatrix/DenseMatrix.hpp>
+
+using DenseMatrix = GridKit::LinearAlgebra::DenseMatrix<double, size_t>;
+
+/**
+ * @brief Class providing methods to evaluate a vector-valued residual and its Jacobian.
+ * This is used to test automatic differentiation.
+ */
+class VectorModel 
+{
+private:
+    std::vector<double> x_, f_;
+    DenseMatrix df_dx_;
+    inline double square_scalar(double);
+    void square(std::vector<double>&, std::vector<double>&);
+
+public:
+    VectorModel(int);
+    void setVariable(std::vector<double>);
+    void evalResidual();
+    void evalJacobian();
+    std::vector<double>& getVariable();
+    std::vector<double>& getResidual();
+    DenseMatrix& getJacobian();
+    ~VectorModel();
+};
diff --git a/examples/Enzyme/Library/library.cpp b/examples/Enzyme/Library/library.cpp
deleted file mode 100644
index 6ff34da7a..000000000
--- a/examples/Enzyme/Library/library.cpp
+++ /dev/null
@@ -1,3 +0,0 @@
-double square(double x) {
-  return x * x;
-}
diff --git a/examples/Enzyme/Library/library.hpp b/examples/Enzyme/Library/library.hpp
deleted file mode 100644
index b498a295b..000000000
--- a/examples/Enzyme/Library/library.hpp
+++ /dev/null
@@ -1,2 +0,0 @@
-#pragma once
-double square(double x);
diff --git a/examples/Enzyme/Library/main.cpp b/examples/Enzyme/Library/main.cpp
deleted file mode 100644
index e0aef4ecd..000000000
--- a/examples/Enzyme/Library/main.cpp
+++ /dev/null
@@ -1,20 +0,0 @@
-#include <iostream>
-#include <limits>
-#include "library.hpp"
-
-double __enzyme_autodiff(double (*)(double), ...);
-
-int main() {
-  int fail = 0;
-  double sq = square(5.0);
-  double dsq = __enzyme_autodiff(square, 5.0);
-
-  std::cout << "x = 5, x^2 = " << sq << ", d(x^2)/dx = " << dsq << "\n"; 
-  if (std::abs(dsq - 10.0) > std::numeric_limits<double>::epsilon())
-  {
-    fail++;
-    std::cout << "Result incorrect\n";
-  }
-  std::cout << "Status: " << fail << "\n";
-  return fail;
-}
diff --git a/examples/Enzyme/PowerElectronics/CMakeLists.txt b/examples/Enzyme/PowerElectronics/CMakeLists.txt
new file mode 100644
index 000000000..a8e74e881
--- /dev/null
+++ b/examples/Enzyme/PowerElectronics/CMakeLists.txt
@@ -0,0 +1,7 @@
+enzyme_add_executable(
+  NAME EnzymePowerElectronicsCheck 
+  SOURCES main.cpp
+  LINK_LIBRARIES GRIDKIT::DenseMatrix GRIDKIT::power_elec_disgen
+)
+
+add_test(NAME "EnzymePowerElectronicsCheck" COMMAND ${CMAKE_CURRENT_BINARY_DIR}/EnzymePowerElectronicsCheck)
diff --git a/examples/Enzyme/PowerElectronics/main.cpp b/examples/Enzyme/PowerElectronics/main.cpp
new file mode 100644
index 000000000..b15b62e5a
--- /dev/null
+++ b/examples/Enzyme/PowerElectronics/main.cpp
@@ -0,0 +1,195 @@
+#include <iostream>
+#include <limits>
+#include <LinearAlgebra/DenseMatrix/DenseMatrix.hpp>
+#include <Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp>
+#include <Model/PowerElectronics/SystemModelPowerElectronics.hpp>
+
+/**
+ * @brief Standalone example that computes the Jacobian associated with the
+ * residual function of DistributedGenerator. We compare the Jacobian obtained
+ * by automatic differentiation via Enzyme to the analytical Jacobian
+ * implemented witin GridKit.
+ *
+ * TODO: Move automatic differentiation inside GridKit and convert this into a unit test.
+ */
+
+using DenseMatrix = GridKit::LinearAlgebra::DenseMatrix<double, size_t>;
+using SparseMatrix = COO_Matrix<double, size_t>;
+using DG = GridKit::DistributedGenerator<double, size_t>;
+using DGParameters = GridKit::DistributedGeneratorParameters<double, size_t>;
+
+int enzyme_dupnoneed;
+int enzyme_dup;
+int enzyme_const;
+void __enzyme_fwddiff(void*, int, std::vector<double>, std::vector<double>,
+                             int, std::vector<double>, std::vector<double>*);
+
+// Copy from DistributedGenerator<ScalarT, IdxT>::evaluateResidual
+// Need to find a way to differentiate the member function directly
+void evaluateResidual(std::vector<double> y_, std::vector<double> f_) 
+{
+    constexpr double wb_ = 2.0*M_PI*50.0;
+    constexpr double wc_ = 31.41;
+    constexpr double mp_ = 9.4e-5;
+    constexpr double Vn_ = 380.0;
+    constexpr double nq_ = 1.3e-3;
+    constexpr double F_ = 0.75;
+    constexpr double Kiv_ = 420.0;
+    constexpr double Kpv_ = 0.1;
+    constexpr double Kic_ = 2.0e4;
+    constexpr double Kpc_ = 15.0;
+    constexpr double Cf_ = 5.0e-5;
+    constexpr double rLf_ = 0.1;
+    constexpr double Lf_ = 1.35e-3;
+    constexpr double rLc_ = 0.03;
+    constexpr double Lc_ = 0.35e-3;
+
+    constexpr bool ref_frame_ = true;
+
+    std::vector<double> yp_(0);
+
+    double omega = wb_ - mp_ * y_[4];
+    if (ref_frame_)
+    {
+        f_[0] = omega - y_[0];
+    }
+    else
+    {
+        f_[0] = 0.0;
+    }
+    
+    f_[1] = cos(y_[3]) * y_[14] - sin(y_[3]) * y_[15];
+    f_[2] = sin(y_[3]) * y_[14] + cos(y_[3]) * y_[15];
+
+    double vbd_in = cos(y_[3]) * y_[1] + sin(y_[3]) * y_[2];
+    double vbq_in = -sin(y_[3]) * y_[1] + cos(y_[3]) * y_[2];
+
+    f_[3] = -yp_[3] + omega - y_[0];
+    f_[4] = -yp_[4] + wc_ * (y_[12] * y_[14] + y_[13] * y_[15] - y_[4]);
+    f_[5] = -yp_[5] + wc_ * (-y_[12] * y_[15] + y_[13] * y_[14] - y_[5]);
+
+    double vod_star = Vn_ - nq_ * y_[5];
+    double voq_star = 0.0;
+
+    f_[6] = -yp_[6] + vod_star - y_[12];
+    f_[7] = -yp_[7] + voq_star - y_[13];
+
+    double ild_star = F_ * y_[14] - wb_ * Cf_ * y_[13] + Kpv_ * (vod_star - y_[12]) + Kiv_ * y_[6];
+    double ilq_star = F_ * y_[15] + wb_ * Cf_ * y_[12] + Kpv_ * (voq_star - y_[13]) + Kiv_ * y_[7];
+
+    f_[8] = -yp_[8] + ild_star - y_[10];
+    f_[9] = -yp_[9] + ilq_star - y_[11];
+
+    double vid_star = -wb_ * Lf_ * y_[11] + Kpc_ * (ild_star - y_[10]) + Kic_ * y_[8];
+    double viq_star = wb_ * Lf_ * y_[10] + Kpc_ * (ilq_star - y_[11]) + Kic_ * y_[9]; 
+    
+    f_[10] = -yp_[10] - (rLf_ / Lf_) * y_[10] + omega * y_[11] +  (vid_star - y_[12]) / Lf_;
+    f_[11] = -yp_[11] - (rLf_ / Lf_) * y_[11] - omega * y_[10] +  (viq_star - y_[13]) / Lf_;
+
+    f_[12] = -yp_[12] + omega * y_[13] + (y_[10] - y_[14]) / Cf_;
+    f_[13] = -yp_[13] - omega * y_[12] + (y_[11] - y_[15]) / Cf_;
+
+    f_[14] = -yp_[14] - (rLc_ / Lc_) * y_[14] + omega * y_[15] + (y_[12] - vbd_in) / Lc_;
+    f_[15] = -yp_[15] - (rLc_ / Lc_) * y_[15] - omega * y_[14] + (y_[13] - vbq_in) / Lc_;
+}
+
+// Function that computes the Jacobian via automatic differentiation
+template <typename T>
+void EnzymeModelJacobian(T* model, DenseMatrix& jac) 
+{
+    int 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)
+    {
+        // Elementary vector for Jacobian-vector product
+        for (int idx = 0; idx < N; ++idx)
+        {
+            y[idx] = (model->y())[idx];
+            res[idx] = (model->getResidual())[idx];
+            v[idx] = 0.0;
+        }
+        v[idy] = 1.0;
+  
+        // Autodiff
+        __enzyme_fwddiff((void*)evaluateResidual, 
+                         enzyme_dup, y, v,
+                         enzyme_dupnoneed, res, &d_res);
+  
+        // Store result
+        for (int idx = 0; idx < N; ++idx)
+        {
+            jac.setValue(idx, idy, d_res[idx]);
+        }
+    }
+}
+
+int main() 
+{
+    // Model
+    DGParameters parms;
+    parms.wb_ = 2.0*M_PI*50.0;
+    parms.wc_ = 31.41;
+    parms.mp_ = 9.4e-5;
+    parms.Vn_ = 380.0;
+    parms.nq_ = 1.3e-3;
+    parms.F_ = 0.75;
+    parms.Kiv_ = 420.0;
+    parms.Kpv_ = 0.1;
+    parms.Kic_ = 2.0e4;
+    parms.Kpc_ = 15.0;
+    parms.Cf_ = 5.0e-5;
+    parms.rLf_ = 0.1;
+    parms.Lf_ = 1.35e-3;
+    parms.rLc_ = 0.03;
+    parms.Lc_ = 0.35e-3;
+    DG *dg = new DG(0, parms, true);
+    dg->allocate();
+    dg->initialize();
+    dg->updateTime(0.0, 0.0);
+
+    // Residual evaluation and reference Jacobian
+    dg->evaluateResidual();
+    dg->evaluateJacobian();
+    std::vector<double> y = dg->y();
+    std::vector<double> yp = dg->yp();
+    std::vector<double> res = dg->getResidual();
+    SparseMatrix jac_ref = dg->getJacobian();
+    DenseMatrix jac_ref_dense(dg->size(), dg->size());
+    jac_ref_dense.setValues(jac_ref);
+  
+    // Enzyme Jacobian
+    DenseMatrix jac_autodiff(dg->size(), dg->size());
+    EnzymeModelJacobian<DG>(dg, jac_autodiff);
+  
+    // Check
+    int fail = 0;
+    bool verbose = true;
+    for (int idy = 0; idy < dg->size(); ++idy)
+    {
+        for (int 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())
+            {
+                fail++;
+                if (verbose)
+                {
+                    std::cout << "Result incorrect at line = " << idy << ", column = " << idx << "\n";
+                }
+            }
+        }
+    }
+    if (verbose)
+    {
+        jac_autodiff.printMatrix("Autodiff Jacobian");
+        jac_ref_dense.printMatrix("Reference Jacobian");
+    }
+    std::cout << "Status: " << fail << "\n";
+
+    // Cleanup
+    delete dg;
+
+    return fail;
+}
diff --git a/examples/Enzyme/Standalone/CMakeLists.txt b/examples/Enzyme/Standalone/CMakeLists.txt
index 919d755e0..fe41dc744 100644
--- a/examples/Enzyme/Standalone/CMakeLists.txt
+++ b/examples/Enzyme/Standalone/CMakeLists.txt
@@ -1,6 +1,13 @@
 enzyme_add_executable(
-  NAME EnzymeStandaloneCheck
-  SOURCES main.cpp
-  )
+  NAME EnzymeStandaloneScalarCheck
+  SOURCES EnzymeScalar.cpp
+)
 
-add_test(NAME "EnzymeStandaloneCheck" COMMAND ${CMAKE_CURRENT_BINARY_DIR}/EnzymeStandaloneCheck)
+enzyme_add_executable(
+  NAME EnzymeStandaloneVectorCheck
+  SOURCES EnzymeVector.cpp
+  LINK_LIBRARIES GRIDKIT::DenseMatrix
+)
+
+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
new file mode 100644
index 000000000..99333b1a8
--- /dev/null
+++ b/examples/Enzyme/Standalone/EnzymeScalar.cpp
@@ -0,0 +1,36 @@
+#include <iostream>
+#include <limits>
+
+/**
+ * @brief Standalone example that computes the derivative of a scalar function
+ * by automatic differentiation via Enzyme.
+ *
+ * TODO: Convert this into a unit test.
+ */
+
+double square(double x) 
+{
+    return x * x;
+}
+
+double __enzyme_autodiff(double(*)(double), ...);
+double dsquare(double x) 
+{
+    return __enzyme_autodiff(square, x);
+}
+
+int main()
+{
+    int fail = 0;
+    double var = 5.0;
+    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())
+    {
+        fail++;
+        std::cout << "Result incorrect\n";
+    }
+    std::cout << "Status: " << fail << "\n";
+    return fail;
+}
diff --git a/examples/Enzyme/Standalone/EnzymeVector.cpp b/examples/Enzyme/Standalone/EnzymeVector.cpp
new file mode 100644
index 000000000..ef6d60c16
--- /dev/null
+++ b/examples/Enzyme/Standalone/EnzymeVector.cpp
@@ -0,0 +1,130 @@
+#include <iostream>
+#include <limits>
+#include <vector>
+#include <LinearAlgebra/DenseMatrix/DenseMatrix.hpp>
+
+/**
+ * @brief Standalone example that computes the Jacobian of a vector-valued function
+ * by automatic differentiation via Enzyme.
+ *
+ * TODO: Convert this into a unit test.
+ */
+
+using DenseMatrix = GridKit::LinearAlgebra::DenseMatrix<double, size_t>;
+int enzyme_dupnoneed;
+int enzyme_dup;
+int enzyme_const;
+void __enzyme_fwddiff(void*, int, std::vector<double>, std::vector<double>, 
+                             int, std::vector<double>, std::vector<double>*);
+
+inline
+double square_scalar(double x) 
+{
+    return x * x;
+}
+
+inline
+double dsquare_ref_scalar(double x) 
+{
+    return 2.0 * x;
+}
+
+// Vector-valued function to differentiate
+void square(std::vector<double> x, std::vector<double>& y) 
+{
+    for (int idx = 0; idx < x.size(); ++idx)
+    {
+        y[idx] = square_scalar(x[idx]);
+    }
+}
+
+// Reference Jacobian
+void dsquare_ref(std::vector<double> x, std::vector<double> y, DenseMatrix& dy) 
+{
+    for (int idy = 0; idy < y.size(); ++idy)
+    {
+        for (int idx = 0; idx < x.size(); ++idx)
+        {
+            if (idx == idy) 
+                dy.setValue(idx, idy, dsquare_ref_scalar(x[idx]));
+        }
+    }
+}
+
+// Function that computes the Jacobian via automatic differentiation
+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)
+    {
+        // Elementary vector for Jacobian-vector product
+        for (int idx = 0; idx < x.size(); ++idx)
+        {
+            v[idx] = 0.0;
+        }
+        v[idy] = 1.0;
+  
+        // Autodiff
+        __enzyme_fwddiff((void*)square, enzyme_dup, x, v,
+                                        enzyme_dupnoneed, y, &d_y);
+  
+        // Store result
+        for (int idx = 0; idx < x.size(); ++idx)
+        {
+            dy.setValue(idx, idy, d_y[idx]);
+        }
+    }
+}
+
+int main()
+{
+    // Vector and matrix declarations
+    constexpr int 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)
+    {
+        x[idx] = rand();
+    }
+  
+    // Function evaluation
+    square(x, sq);
+  
+    // Reference Jacobian
+    dsquare_ref(x, sq, dsq_ref);
+  
+    // Enzyme Jacobian
+    dsquare(x, sq, dsq);
+  
+    // Check
+    int fail = 0;
+    bool verbose = true;
+    for (int idy = 0; idy < sq.size(); ++idy)
+    {
+        for (int idx = 0; idx < x.size(); ++idx)
+        {
+            if (std::abs(dsq.getValue(idx, idy) - dsq_ref.getValue(idx, idy)) > std::numeric_limits<double>::epsilon())
+            {
+                fail++;
+                if (verbose)
+                {
+                    std::cout << "Result incorrect at line = " << idy << ", column = " << idx << "\n";
+                    std::cout << "x = " << x[idx] << ", x^2 = " << sq[idx] << ", d(x^2)/dx = " << dsq.getValue(idx, idy) << "\n"; 
+                }
+            }
+        }
+    }
+    if (verbose)
+    {
+        dsq.printMatrix("Autodiff Jacobian");
+        dsq_ref.printMatrix("Reference Jacobian");
+    }
+    std::cout << "Status: " << fail << "\n";
+    return fail;
+}
diff --git a/examples/Enzyme/Standalone/main.cpp b/examples/Enzyme/Standalone/main.cpp
deleted file mode 100644
index c9b81ef05..000000000
--- a/examples/Enzyme/Standalone/main.cpp
+++ /dev/null
@@ -1,27 +0,0 @@
-#include <iostream>
-#include <limits>
-
-double square(double x) {
-  return x * x;
-}
-
-double __enzyme_autodiff(double(*)(double), ...);
-double dsquare(double x) {
-  return __enzyme_autodiff(square, x);
-}
-
-double sq  = square(5.0);
-double dsq = dsquare(5.0);
-
-int main()
-{
-  int fail = 0;
-  std::cout << "x = 5, x^2 = " << sq << ", d(x^2)/dx = " << dsq << "\n"; 
-  if (std::abs(dsq - 10.0) > std::numeric_limits<double>::epsilon())
-  {
-    fail++;
-    std::cout << "Result incorrect\n";
-  }
-  std::cout << "Status: " << fail << "\n";
-  return fail;
-}
diff --git a/examples/LinearAlgebra/CMakeLists.txt b/examples/LinearAlgebra/CMakeLists.txt
new file mode 100644
index 000000000..553b0587e
--- /dev/null
+++ b/examples/LinearAlgebra/CMakeLists.txt
@@ -0,0 +1,3 @@
+
+add_subdirectory(SparseTest)
+add_subdirectory(DenseTest)
diff --git a/examples/LinearAlgebra/DenseTest/CMakeLists.txt b/examples/LinearAlgebra/DenseTest/CMakeLists.txt
new file mode 100644
index 000000000..5d3bf7b5c
--- /dev/null
+++ b/examples/LinearAlgebra/DenseTest/CMakeLists.txt
@@ -0,0 +1,6 @@
+
+add_executable(densemattest DenseTest.cpp)
+target_link_libraries(densemattest GRIDKIT::DenseMatrix)
+                                
+add_test(NAME DenseMatrixTest COMMAND $<TARGET_FILE:densemattest>)
+install(TARGETS densemattest RUNTIME DESTINATION bin)
diff --git a/examples/LinearAlgebra/DenseTest/DenseTest.cpp b/examples/LinearAlgebra/DenseTest/DenseTest.cpp
new file mode 100644
index 000000000..8b5e27bd3
--- /dev/null
+++ b/examples/LinearAlgebra/DenseTest/DenseTest.cpp
@@ -0,0 +1,26 @@
+
+#include <LinearAlgebra/DenseMatrix/DenseMatrix.hpp>
+
+int main()
+{
+    int fail = 0;
+
+    size_t m = 4;
+    size_t n = 4;
+    GridKit::LinearAlgebra::DenseMatrix<double, size_t> A = 
+      GridKit::LinearAlgebra::DenseMatrix<double, size_t>(m, n);
+
+    double val = 0.0;
+    for (size_t j = 0; j < n; ++j)
+    {
+        for (size_t i = 0; i < m; ++i)
+        {
+            A.setValue(i, j, val);
+            val += 1.0;
+        }
+    }
+    A.printMatrix("Dense matrix test output");
+    (A.getValuesCOO())->printMatrix();
+
+    return fail;
+}
diff --git a/examples/SparseTest/CMakeLists.txt b/examples/LinearAlgebra/SparseTest/CMakeLists.txt
similarity index 100%
rename from examples/SparseTest/CMakeLists.txt
rename to examples/LinearAlgebra/SparseTest/CMakeLists.txt
diff --git a/examples/SparseTest/SparseTest.cpp b/examples/LinearAlgebra/SparseTest/SparseTest.cpp
similarity index 97%
rename from examples/SparseTest/SparseTest.cpp
rename to examples/LinearAlgebra/SparseTest/SparseTest.cpp
index 32c075c45..56a70ea15 100644
--- a/examples/SparseTest/SparseTest.cpp
+++ b/examples/LinearAlgebra/SparseTest/SparseTest.cpp
@@ -8,7 +8,7 @@
 #include <vector>
 #include <tuple>
 
-#include <LinearAlgebra/COO_Matrix.hpp>
+#include <LinearAlgebra/SparseMatrix/COO_Matrix.hpp>
 
 int main(int argc, char const *argv[])
 {
diff --git a/examples/Microgrid/CMakeLists.txt b/examples/Microgrid/CMakeLists.txt
index 77c59d753..9d627fcac 100644
--- a/examples/Microgrid/CMakeLists.txt
+++ b/examples/Microgrid/CMakeLists.txt
@@ -3,11 +3,11 @@
 
 
 add_executable(microgrid Microgrid.cpp)
-target_link_libraries(microgrid GRIDKIT::powerelec_disgen
-                                GRIDKIT::powerelec_mircoline
-                                GRIDKIT::powerelec_microload
+target_link_libraries(microgrid GRIDKIT::power_elec_disgen
+                                GRIDKIT::power_elec_microline
+                                GRIDKIT::power_elec_microload
                                 GRIDKIT::solvers_dyn
-                                GRIDKIT::powerelec_mircobusdq)
+                                GRIDKIT::power_elec_microbusdq)
                                 
 add_test(NAME Microgrid    COMMAND $<TARGET_FILE:microgrid>)
 install(TARGETS microgrid RUNTIME DESTINATION bin)
diff --git a/examples/RLCircuit/CMakeLists.txt b/examples/RLCircuit/CMakeLists.txt
index e7e052603..56599c521 100644
--- a/examples/RLCircuit/CMakeLists.txt
+++ b/examples/RLCircuit/CMakeLists.txt
@@ -3,10 +3,10 @@
 
 
 add_executable(rlcircuit RLCircuit.cpp)
-target_link_libraries(rlcircuit GRIDKIT::powerelec_capacitor
-                                GRIDKIT::powerelec_inductor
-                                GRIDKIT::powerelec_resistor
-                                GRIDKIT::powerelec_voltagesource
+target_link_libraries(rlcircuit GRIDKIT::power_elec_capacitor
+                                GRIDKIT::power_elec_inductor
+                                GRIDKIT::power_elec_resistor
+                                GRIDKIT::power_elec_voltagesource
                                 GRIDKIT::solvers_dyn)
                                 
 add_test(NAME RLCircuit    COMMAND $<TARGET_FILE:rlcircuit>)
diff --git a/examples/ScaleMicrogrid/CMakeLists.txt b/examples/ScaleMicrogrid/CMakeLists.txt
index d1822c587..1f0fe0cbb 100644
--- a/examples/ScaleMicrogrid/CMakeLists.txt
+++ b/examples/ScaleMicrogrid/CMakeLists.txt
@@ -3,11 +3,11 @@
 
 
 add_executable(scalemicrogrid ScaleMicrogrid.cpp)
-target_link_libraries(scalemicrogrid GRIDKIT::powerelec_disgen
-                                     GRIDKIT::powerelec_mircoline
-                                     GRIDKIT::powerelec_microload
+target_link_libraries(scalemicrogrid GRIDKIT::power_elec_disgen
+                                     GRIDKIT::power_elec_microline
+                                     GRIDKIT::power_elec_microload
                                      GRIDKIT::solvers_dyn
-                                     GRIDKIT::powerelec_mircobusdq)
+                                     GRIDKIT::power_elec_microbusdq)
 								
 add_test(NAME ScaleMicrogrid    COMMAND $<TARGET_FILE:scalemicrogrid>)
 install(TARGETS scalemicrogrid RUNTIME DESTINATION bin)
diff --git a/src/LinearAlgebra/CMakeLists.txt b/src/LinearAlgebra/CMakeLists.txt
index b51d669ec..d290617be 100644
--- a/src/LinearAlgebra/CMakeLists.txt
+++ b/src/LinearAlgebra/CMakeLists.txt
@@ -1,7 +1,3 @@
 
-add_library(SparseMatrix INTERFACE)
-include_directories(SparseMatrix INTERFACE ${CMAKE_CURRENT_LIST_DIR})
-
-add_library(GRIDKIT::SparseMatrix ALIAS SparseMatrix)
-
-
+add_subdirectory(SparseMatrix)
+add_subdirectory(DenseMatrix)
diff --git a/src/LinearAlgebra/DenseMatrix/CMakeLists.txt b/src/LinearAlgebra/DenseMatrix/CMakeLists.txt
new file mode 100644
index 000000000..37eb294c3
--- /dev/null
+++ b/src/LinearAlgebra/DenseMatrix/CMakeLists.txt
@@ -0,0 +1,6 @@
+
+add_library(DenseMatrix INTERFACE)
+include_directories(DenseMatrix INTERFACE ${CMAKE_CURRENT_LIST_DIR})
+
+add_library(GRIDKIT::DenseMatrix ALIAS DenseMatrix)
+
diff --git a/src/LinearAlgebra/DenseMatrix/DenseMatrix.hpp b/src/LinearAlgebra/DenseMatrix/DenseMatrix.hpp
new file mode 100644
index 000000000..839156f92
--- /dev/null
+++ b/src/LinearAlgebra/DenseMatrix/DenseMatrix.hpp
@@ -0,0 +1,230 @@
+#pragma once
+
+#include <iostream>
+#include <limits>
+#include <vector>
+#include <assert.h>
+#include <LinearAlgebra/SparseMatrix/COO_Matrix.hpp>
+
+/**
+ * @brief Class to provide dense matrices. 
+ *
+ * This is intended for small matrices that store model Jacobians to be subsequently copied 
+ * into large sparse matrices.
+ */
+namespace GridKit
+{
+namespace LinearAlgebra
+{
+template <class ScalarT, typename IdxT>
+class DenseMatrix
+{
+private:
+    IdxT rows_size_;
+    IdxT columns_size_;
+    std::vector<ScalarT> values_;
+    COO_Matrix<ScalarT, IdxT> values_COO_;
+    bool values_changed_ = false;
+    bool sparsified_ = false;
+public:
+    // Constructors and destructors
+    DenseMatrix(const IdxT rows_size, const IdxT columns_size);
+    ~DenseMatrix();
+  
+    // Getters and setters
+    ScalarT getValue(const IdxT i, const IdxT j) const;
+    void setValue(const IdxT i, const IdxT j, const ScalarT value);
+    void setValues(COO_Matrix<ScalarT, IdxT> values_COO);
+    std::vector<ScalarT>* getValues();
+    COO_Matrix<ScalarT, IdxT>* getValuesCOO();
+
+    // Utilities
+    void toCOO();
+    void printMatrix(std::string name="");
+
+    // Purposefully not defining BLAS operations. This class should not be used
+    // for compute.
+};
+
+/**
+ * @brief DenseMatrix constructor
+ *
+ * @tparam ScalarT 
+ * @tparam IdxT 
+ *
+ * @param[in] IdxT - rows_size 
+ * @param[in] IdxT - columns_size 
+ */
+template <class ScalarT, typename IdxT>
+DenseMatrix<ScalarT, IdxT>::DenseMatrix(const IdxT rows_size, const IdxT columns_size) :
+    rows_size_(rows_size),
+    columns_size_(columns_size),
+    values_(rows_size*columns_size, 0),
+    values_COO_(rows_size, columns_size)
+{
+
+}
+
+/**
+ * @brief DenseMatrix single value getter
+ *
+ * @tparam ScalarT 
+ * @tparam IdxT 
+ *
+ * @param[in] IdxT - i row index
+ * @param[in] IdxT - j column index 
+ * @return ScalarT - value
+ */
+template <class ScalarT, typename IdxT>
+inline ScalarT DenseMatrix<ScalarT, IdxT>::getValue(const IdxT i, const IdxT j) const
+{
+    assert(i < this->columns_size_);
+    assert(j < this->rows_size_);
+    return this->values_[j*rows_size_+i];
+}
+
+/**
+ * @brief DenseMatrix single value setter
+ *
+ * @tparam ScalarT 
+ * @tparam IdxT 
+ *
+ * @param[in] IdxT - i row index
+ * @param[in] IdxT - j column index 
+ * @param[in] ScalarT - value 
+ */
+template <class ScalarT, typename IdxT>
+inline void DenseMatrix<ScalarT, IdxT>::setValue(const IdxT i, const IdxT j, const ScalarT value)
+{
+    assert(i < this->columns_size_);
+    assert(j < this->rows_size_);
+    this->values_[j*rows_size_+i] = value;
+    values_changed_ = true;
+}
+
+/**
+ * @brief DenseMatrix value setter from COO
+ *
+ * @tparam ScalarT 
+ * @tparam IdxT 
+ *
+ * @param[in] COO_Matrix<ScalarT, IdxT> - values_COO
+ */
+template <class ScalarT, typename IdxT>
+inline void DenseMatrix<ScalarT, IdxT>::setValues(COO_Matrix<ScalarT, IdxT> values_COO)
+{
+    std::tuple<std::vector<IdxT>&, std::vector<IdxT>&, std::vector<ScalarT>&> entries = values_COO.getEntries();
+    const auto [rcord, ccord, vals] = entries;
+    for (IdxT idx = 0; idx < values_COO.nnz(); ++idx)
+    {
+        this->setValue(rcord[idx], ccord[idx], vals[idx]);
+    }
+}
+
+/**
+ * @brief DenseMatrix getter for all values stored as a vector
+ *
+ * @tparam ScalarT 
+ * @tparam IdxT 
+ *
+ * @return Address of the vector containing matrix values
+ */
+template <class ScalarT, typename IdxT>
+inline std::vector<ScalarT>* DenseMatrix<ScalarT, IdxT>::getValues()
+{
+    return &(this->values_);
+}
+
+/**
+ * @brief DenseMatrix getter for all values stored as a COO sparse matrix
+ *
+ * @tparam ScalarT 
+ * @tparam IdxT 
+ *
+ * @return Address of the COO matrix containing the sparsified matrix values
+ */
+template <class ScalarT, typename IdxT>
+inline COO_Matrix<ScalarT, IdxT>* DenseMatrix<ScalarT, IdxT>::getValuesCOO()
+{
+    if (!sparsified_ || values_changed_)
+    {
+        this->toCOO();
+    }
+    return &(this->values_COO_);
+}
+
+/**
+ * @brief Dense matrix conversion to COO form
+ *
+ * @tparam ScalarT 
+ * @tparam IdxT 
+ */
+template <class ScalarT, typename IdxT>
+inline void DenseMatrix<ScalarT, IdxT>::toCOO()
+{
+    if (!sparsified_ || values_changed_)
+    {
+        IdxT nnz = 0;
+        std::vector<IdxT> rcord;
+        std::vector<IdxT> ccord;
+        std::vector<ScalarT> vals;
+        for (IdxT j = 0; j < this->columns_size_; ++j)
+        {
+            for (IdxT i = 0; i < this->rows_size_; ++i)
+            {
+                ScalarT value = this->values_[j*rows_size_+i];
+                if (std::abs(value) > std::numeric_limits<double>::epsilon())
+                {
+                    nnz++;
+                    rcord.push_back(i); 
+                    ccord.push_back(j); 
+                    vals.push_back(value); 
+                }
+            }
+        }
+        values_COO_.setValues(rcord, ccord, vals);
+        sparsified_ = true;
+        values_changed_ = false;
+    }
+}
+
+/**
+ * @brief Print matrix
+ *
+ * @tparam ScalarT 
+ * @tparam IdxT 
+ *
+ * @param[in] name to identify the specific matrix printed
+ */
+template <class ScalarT, typename IdxT>
+inline void DenseMatrix<ScalarT, IdxT>::printMatrix(std::string name)
+{
+    std::cout << "Dense matrix: " << name << "\n";
+    for (IdxT i = 0; i < this->rows_size_; ++i)
+    {
+        for (IdxT j = 0; j < this->columns_size_; ++j)
+        {
+          std::cout << this->values_[j*rows_size_+i] << " ";
+        }
+        std::cout << "\n";
+    }
+}
+
+/**
+ * @brief DenseMatrix destructor
+ *
+ * @tparam ScalarT 
+ * @tparam IdxT 
+ */
+template <class ScalarT, typename IdxT>
+DenseMatrix<ScalarT, IdxT>::~DenseMatrix()
+{
+    
+}
+
+// Available template instantiations
+template class DenseMatrix<double, size_t>;
+
+} // LinearAlgebra
+} // GridKit
+
diff --git a/src/LinearAlgebra/SparseMatrix/CMakeLists.txt b/src/LinearAlgebra/SparseMatrix/CMakeLists.txt
new file mode 100644
index 000000000..ddb185489
--- /dev/null
+++ b/src/LinearAlgebra/SparseMatrix/CMakeLists.txt
@@ -0,0 +1,6 @@
+
+add_library(SparseMatrix INTERFACE)
+include_directories(SparseMatrix INTERFACE ${CMAKE_CURRENT_LIST_DIR})
+
+add_library(GRIDKIT::SparseMatrix ALIAS SparseMatrix)
+
diff --git a/src/LinearAlgebra/COO_Matrix.hpp b/src/LinearAlgebra/SparseMatrix/COO_Matrix.hpp
similarity index 98%
rename from src/LinearAlgebra/COO_Matrix.hpp
rename to src/LinearAlgebra/SparseMatrix/COO_Matrix.hpp
index b84e239e1..ea72e41a5 100644
--- a/src/LinearAlgebra/COO_Matrix.hpp
+++ b/src/LinearAlgebra/SparseMatrix/COO_Matrix.hpp
@@ -69,7 +69,7 @@ class COO_Matrix
 
     std::tuple<IdxT, IdxT> getDimensions();
 
-    void printMatrix();
+    void printMatrix(std::string name="");
 
     
     static void sortSparseCOO(std::vector<IdxT> &rows, std::vector<IdxT> &columns, std::vector<ScalarT> &values);
@@ -443,7 +443,6 @@ inline ScalarT COO_Matrix<ScalarT, IdxT>::frobNorm()
  * @brief Permutate the matrix to a different one. Only changes the coordinates
  * 
  * @tparam ScalarT 
- * @tparam IdxT 
  * @param[in] row_perm 
  * @param[out] col_perm 
  * 
@@ -589,7 +588,7 @@ inline bool COO_Matrix<ScalarT, IdxT>::isSorted()
 template <class ScalarT, typename IdxT>
 inline IdxT COO_Matrix<ScalarT, IdxT>::nnz()
 {
-    return static_cast<IdxT>(this->values_.size);
+    return static_cast<IdxT>(this->values_.size());
 }
 
 template <class ScalarT, typename IdxT>
@@ -603,16 +602,18 @@ inline std::tuple<IdxT, IdxT> COO_Matrix<ScalarT, IdxT>::getDimensions()
  * 
  * @tparam ScalarT 
  * @tparam IdxT 
+ * 
+ * @param[in] name to identify the specific matrix printed
  */
 template <class ScalarT, typename IdxT>
-inline void COO_Matrix<ScalarT, IdxT>::printMatrix()
+inline void COO_Matrix<ScalarT, IdxT>::printMatrix(std::string name)
 {
     if (this->sorted_ == false)
     {
         this->sortSparse();
     }
     
-    std::cout << "Sparse COO Matrix\n";
+    std::cout << "Sparse COO Matrix: " << name << "\n";
     std::cout << "(x , y, value)\n";
     for (size_t i = 0; i < this->values_.size(); i++)
     {
@@ -879,4 +880,4 @@ COO_Matrix<ScalarT, IdxT>::~COO_Matrix()
     
 }
 
-#endif
\ No newline at end of file
+#endif
diff --git a/src/Model/Evaluator.hpp b/src/Model/Evaluator.hpp
index 669834dcf..7a374d975 100644
--- a/src/Model/Evaluator.hpp
+++ b/src/Model/Evaluator.hpp
@@ -2,7 +2,7 @@
 
 #include <vector>
 #include <ScalarTraits.hpp>
-#include <LinearAlgebra/COO_Matrix.hpp>
+#include <LinearAlgebra/SparseMatrix/COO_Matrix.hpp>
 
 namespace GridKit
 {
diff --git a/src/Model/PowerElectronics/Capacitor/CMakeLists.txt b/src/Model/PowerElectronics/Capacitor/CMakeLists.txt
index 62efc1795..2e05fad26 100644
--- a/src/Model/PowerElectronics/Capacitor/CMakeLists.txt
+++ b/src/Model/PowerElectronics/Capacitor/CMakeLists.txt
@@ -1,8 +1,8 @@
 
 
 
-gridkit_add_library(powerelec_capacitor
+gridkit_add_library(power_elec_capacitor
   SOURCES
     Capacitor.cpp
   OUTPUT_NAME
-    gridkit_powerelec_capacitor)
\ No newline at end of file
+    gridkit_power_elec_capacitor)
\ No newline at end of file
diff --git a/src/Model/PowerElectronics/DistributedGenerator/CMakeLists.txt b/src/Model/PowerElectronics/DistributedGenerator/CMakeLists.txt
index 22652f53b..3d4b805cb 100644
--- a/src/Model/PowerElectronics/DistributedGenerator/CMakeLists.txt
+++ b/src/Model/PowerElectronics/DistributedGenerator/CMakeLists.txt
@@ -1,8 +1,8 @@
 
 
 
-gridkit_add_library(powerelec_disgen
+gridkit_add_library(power_elec_disgen
   SOURCES
     DistributedGenerator.cpp
   OUTPUT_NAME
-    gridkit_powerelec_disgen)
\ No newline at end of file
+    gridkit_power_elec_disgen)
\ No newline at end of file
diff --git a/src/Model/PowerElectronics/InductionMotor/CMakeLists.txt b/src/Model/PowerElectronics/InductionMotor/CMakeLists.txt
index 95ef42436..ce94553a9 100644
--- a/src/Model/PowerElectronics/InductionMotor/CMakeLists.txt
+++ b/src/Model/PowerElectronics/InductionMotor/CMakeLists.txt
@@ -1,8 +1,8 @@
 
 
 
-gridkit_add_library(powerelec_inductionmotor
+gridkit_add_library(power_elec_inductionmotor
   SOURCES
     InductionMotor.cpp
   OUTPUT_NAME
-    gridkit_powerelec_inductionmotor)
\ No newline at end of file
+    gridkit_power_elec_inductionmotor)
\ No newline at end of file
diff --git a/src/Model/PowerElectronics/Inductor/CMakeLists.txt b/src/Model/PowerElectronics/Inductor/CMakeLists.txt
index 4c620030f..3570a451f 100644
--- a/src/Model/PowerElectronics/Inductor/CMakeLists.txt
+++ b/src/Model/PowerElectronics/Inductor/CMakeLists.txt
@@ -1,8 +1,8 @@
 
 
 
-gridkit_add_library(powerelec_inductor
+gridkit_add_library(power_elec_inductor
   SOURCES
     Inductor.cpp
   OUTPUT_NAME
-    gridkit_powerelec_inductor)
\ No newline at end of file
+    gridkit_power_elec_inductor)
\ No newline at end of file
diff --git a/src/Model/PowerElectronics/LinearTransformer/CMakeLists.txt b/src/Model/PowerElectronics/LinearTransformer/CMakeLists.txt
index eaf95c043..da1e7efba 100644
--- a/src/Model/PowerElectronics/LinearTransformer/CMakeLists.txt
+++ b/src/Model/PowerElectronics/LinearTransformer/CMakeLists.txt
@@ -1,8 +1,8 @@
 
 
 
-gridkit_add_library(powerelec_lineartrasnformer
+gridkit_add_library(power_elec_lineartrasnformer
   SOURCES
     LinearTransformer.cpp
   OUTPUT_NAME
-    gridkit_powerelec_lineartrasnformer)
\ No newline at end of file
+    gridkit_power_elec_lineartrasnformer)
\ No newline at end of file
diff --git a/src/Model/PowerElectronics/MicrogridBusDQ/CMakeLists.txt b/src/Model/PowerElectronics/MicrogridBusDQ/CMakeLists.txt
index 1bfb8dc3d..21efbb7d4 100644
--- a/src/Model/PowerElectronics/MicrogridBusDQ/CMakeLists.txt
+++ b/src/Model/PowerElectronics/MicrogridBusDQ/CMakeLists.txt
@@ -1,8 +1,8 @@
 
 
 
-gridkit_add_library(powerelec_mircobusdq
+gridkit_add_library(power_elec_microbusdq
   SOURCES
     MicrogridBusDQ.cpp
   OUTPUT_NAME
-    gridkit_powerelec_mircobusdq)
\ No newline at end of file
+    gridkit_power_elec_microbusdq)
\ No newline at end of file
diff --git a/src/Model/PowerElectronics/MicrogridLine/CMakeLists.txt b/src/Model/PowerElectronics/MicrogridLine/CMakeLists.txt
index d59cebb9c..3f4224600 100644
--- a/src/Model/PowerElectronics/MicrogridLine/CMakeLists.txt
+++ b/src/Model/PowerElectronics/MicrogridLine/CMakeLists.txt
@@ -1,8 +1,8 @@
 
 
 
-gridkit_add_library(powerelec_mircoline
+gridkit_add_library(power_elec_microline
   SOURCES
     MicrogridLine.cpp
   OUTPUT_NAME
-    gridkit_powerelec_mircoline)
\ No newline at end of file
+    gridkit_power_elec_microline)
\ No newline at end of file
diff --git a/src/Model/PowerElectronics/MicrogridLoad/CMakeLists.txt b/src/Model/PowerElectronics/MicrogridLoad/CMakeLists.txt
index 313d47688..7f9eaf671 100644
--- a/src/Model/PowerElectronics/MicrogridLoad/CMakeLists.txt
+++ b/src/Model/PowerElectronics/MicrogridLoad/CMakeLists.txt
@@ -1,8 +1,8 @@
 
 
 
-gridkit_add_library(powerelec_microload
+gridkit_add_library(power_elec_microload
   SOURCES
     MicrogridLoad.cpp
   OUTPUT_NAME
-    gridkit_powerelec_microload)
\ No newline at end of file
+    gridkit_power_elec_microload)
\ No newline at end of file
diff --git a/src/Model/PowerElectronics/Resistor/CMakeLists.txt b/src/Model/PowerElectronics/Resistor/CMakeLists.txt
index 9386bda83..757597a6a 100644
--- a/src/Model/PowerElectronics/Resistor/CMakeLists.txt
+++ b/src/Model/PowerElectronics/Resistor/CMakeLists.txt
@@ -1,8 +1,8 @@
 
 
 
-gridkit_add_library(powerelec_resistor
+gridkit_add_library(power_elec_resistor
   SOURCES
     Resistor.cpp
   OUTPUT_NAME
-    gridkit_powerelec_resistor)
\ No newline at end of file
+    gridkit_power_elec_resistor)
\ No newline at end of file
diff --git a/src/Model/PowerElectronics/SynchronousMachine/CMakeLists.txt b/src/Model/PowerElectronics/SynchronousMachine/CMakeLists.txt
index bfaeace4b..c42816d96 100644
--- a/src/Model/PowerElectronics/SynchronousMachine/CMakeLists.txt
+++ b/src/Model/PowerElectronics/SynchronousMachine/CMakeLists.txt
@@ -1,8 +1,8 @@
 
 
 
-gridkit_add_library(powerelec_synmachine
+gridkit_add_library(power_elec_synmachine
   SOURCES
     SynchronousMachine.cpp
   OUTPUT_NAME
-    gridkit_powerelec_synmachine)
\ No newline at end of file
+    gridkit_power_elec_synmachine)
\ No newline at end of file
diff --git a/src/Model/PowerElectronics/TransmissionLine/CMakeLists.txt b/src/Model/PowerElectronics/TransmissionLine/CMakeLists.txt
index 6e2cea4f4..c80b3b6de 100644
--- a/src/Model/PowerElectronics/TransmissionLine/CMakeLists.txt
+++ b/src/Model/PowerElectronics/TransmissionLine/CMakeLists.txt
@@ -1,8 +1,8 @@
 
 
 
-gridkit_add_library(powerelec_tranline
+gridkit_add_library(power_elec_tranline
   SOURCES
     TransmissionLine.cpp
   OUTPUT_NAME
-    gridkit_powerelec_tranline)
\ No newline at end of file
+    gridkit_power_elec_tranline)
\ No newline at end of file
diff --git a/src/Model/PowerElectronics/VoltageSource/CMakeLists.txt b/src/Model/PowerElectronics/VoltageSource/CMakeLists.txt
index 7196f4d43..4fc504742 100644
--- a/src/Model/PowerElectronics/VoltageSource/CMakeLists.txt
+++ b/src/Model/PowerElectronics/VoltageSource/CMakeLists.txt
@@ -1,8 +1,8 @@
 
 
 
-gridkit_add_library(powerelec_voltagesource
+gridkit_add_library(power_elec_voltagesource
   SOURCES
   VoltageSource.cpp
   OUTPUT_NAME
-    gridkit_powerelec_voltagesource)
\ No newline at end of file
+    gridkit_power_elec_voltagesource)
\ No newline at end of file