diff --git a/GridKit/LinearAlgebra/Vector/Vector.cpp b/GridKit/LinearAlgebra/Vector/Vector.cpp
index da10f393d..e5aaa77b8 100644
--- a/GridKit/LinearAlgebra/Vector/Vector.cpp
+++ b/GridKit/LinearAlgebra/Vector/Vector.cpp
@@ -12,6 +12,37 @@ namespace GridKit
 
     using out = GridKit::Utilities::Logger;
 
+    namespace detail
+    {
+      [[gnu::cold, gnu::noinline]]
+      void logHostUnsyncFailure()
+      {
+        out::error() << "Vector::getData - host data is stale. Perhaps you need to call syncData?\n";
+      }
+
+      [[gnu::cold, gnu::noinline]]
+      void logDeviceUnsyncFailure()
+      {
+        out::error() << "Vector::getData - host device is stale. Perhaps you need to call syncData?\n";
+      }
+
+      template <typename IdxT>
+      [[gnu::cold, gnu::noinline]]
+      void logBoundsCheckFailure(IdxT j, IdxT k)
+      {
+        out::error() << "Vector::getData - vector index " << j << " out of range, multivector has only " << k << " vectors\n";
+      }
+
+      template <typename IdxT>
+      [[gnu::cold, gnu::noinline]]
+      void logUpdatedBoundsCheckFailure(IdxT j, IdxT k)
+      {
+        out::error() << "Vector::setDataUpdated - vector index " << j
+                     << " out of range, multivector has only " << k
+                     << " vectors\n";
+      }
+    } // namespace detail
+
     /**
      * @brief Single vector constructor.
      *
@@ -22,11 +53,9 @@ namespace GridKit
       : n_capacity_(n),
         k_(1),
         n_size_(n),
-        gpu_updated_(new bool[1]),
-        cpu_updated_(new bool[1])
+        gpu_updated_spec_(false),
+        cpu_updated_spec_(false)
     {
-      gpu_updated_[0] = false;
-      cpu_updated_[0] = false;
     }
 
     /**
@@ -39,12 +68,21 @@ namespace GridKit
     Vector<ScalarT, IdxT>::Vector(IdxT n, IdxT k)
       : n_capacity_(n),
         k_(k),
-        n_size_(n),
-        gpu_updated_(new bool[static_cast<std::size_t>(k)]),
-        cpu_updated_(new bool[static_cast<std::size_t>(k)])
+        n_size_(n)
     {
-      setHostUpdated(false);
-      setDeviceUpdated(false);
+      if (k <= 1)
+      {
+        gpu_updated_spec_ = false;
+        cpu_updated_spec_ = false;
+      }
+      else
+      {
+        gpu_updated_ = new bool[static_cast<std::size_t>(k)];
+        cpu_updated_ = new bool[static_cast<std::size_t>(k)];
+
+        setHostUpdated(false);
+        setDeviceUpdated(false);
+      }
     }
 
     /**
@@ -58,8 +96,12 @@ namespace GridKit
         mem_.deleteOnHost(h_data_);
       if (owns_gpu_data_ && d_data_)
         mem_.deleteOnDevice(d_data_);
-      delete[] gpu_updated_;
-      delete[] cpu_updated_;
+
+      if (k_ > 1)
+      {
+        delete[] gpu_updated_;
+        delete[] cpu_updated_;
+      }
     }
 
     /**
@@ -217,79 +259,6 @@ namespace GridKit
       return 0;
     }
 
-    /**
-     * @brief Set the flag to indicate that the data (HOST or DEVICE) has been
-     * updated.
-     *
-     * Use this function if you update vector elements by accessing the raw data
-     * pointer.
-     *
-     * @param[in] memspace - Memory space (HOST or DEVICE)
-     *
-     * @warning This is an expert level method. Use only if you know what
-     * you are doing.
-     */
-    template <typename ScalarT, typename IdxT>
-    int Vector<ScalarT, IdxT>::setDataUpdated(memory::MemorySpace memspace)
-    {
-      assert(cpu_updated_ && gpu_updated_ && "Update flags not allocated");
-
-      using namespace memory;
-      switch (memspace)
-      {
-      case HOST:
-        setHostUpdated(true);
-        setDeviceUpdated(false);
-        break;
-      case DEVICE:
-        setHostUpdated(false);
-        setDeviceUpdated(true);
-        break;
-      }
-      return 0;
-    }
-
-    /**
-     * @brief Set the flag to indicate that the data (HOST or DEVICE) for
-     * vector `j` in the multivector has been updated.
-     *
-     * Use this function if you update vector elements by accessing the raw data
-     * pointer.
-     *
-     * @param[in] memspace - Memory space (HOST or DEVICE)
-     *
-     * @warning This is an expert level method. Use only if you know what
-     * you are doing.
-     */
-    template <typename ScalarT, typename IdxT>
-    int Vector<ScalarT, IdxT>::setDataUpdated(IdxT j, memory::MemorySpace memspace)
-    {
-      assert(cpu_updated_ && gpu_updated_ && "Update flags not allocated");
-
-      using namespace memory;
-
-      if (k_ <= j)
-      {
-        out::error() << "Vector::setDataUpdated - vector index " << j
-                     << " out of range, multivector has only " << k_
-                     << " vectors\n";
-        return 1;
-      }
-
-      switch (memspace)
-      {
-      case HOST:
-        cpu_updated_[j] = true;
-        gpu_updated_[j] = false;
-        break;
-      case DEVICE:
-        gpu_updated_[j] = true;
-        cpu_updated_[j] = false;
-        break;
-      }
-      return 0;
-    }
-
     /**
      * @brief Copy data from another vector.
      *
@@ -389,171 +358,6 @@ namespace GridKit
       return 0;
     }
 
-    /**
-     * @brief get a pointer to HOST or DEVICE vector data.
-     *
-     * @param[in] memspace  - Memory space of the pointer (HOST or DEVICE)
-     *
-     * @return pointer to the vector data (HOST or DEVICE). In case of multivectors,
-     * vectors are stored column-wise.
-     *
-     * @note This function gives you access to the pointer, not to a copy.
-     * If you change the values using the pointer, the vector values will
-     * change too. Make sure to use setDataUpdated function to set the update
-     * flags correctly after changing the values.
-     */
-    template <typename ScalarT, typename IdxT>
-    ScalarT* Vector<ScalarT, IdxT>::getData(memory::MemorySpace memspace)
-    {
-      using memory::DEVICE;
-      using memory::HOST;
-
-      switch (memspace)
-      {
-      case HOST:
-        if (cpu_updated_[0] == false)
-        {
-          out::error() << "Vector::getData - host data is stale. Perhaps you need to call syncData?\n";
-          return nullptr;
-        }
-        return h_data_;
-      case DEVICE:
-        if (gpu_updated_[0] == false)
-        {
-          out::error() << "Vector::getData - device data is stale. Perhaps you need to call syncData?\n";
-          return nullptr;
-        }
-        return d_data_;
-      default:
-        return nullptr;
-      }
-    }
-
-    /**
-     * @brief get a pointer to HOST or DEVICE vector data.
-     *
-     * @param[in] memspace  - Memory space of the pointer (HOST or DEVICE)
-     *
-     * @return pointer to the vector data (HOST or DEVICE). In case of multivectors,
-     * vectors are stored column-wise.
-     */
-    template <typename ScalarT, typename IdxT>
-    const ScalarT* Vector<ScalarT, IdxT>::getData(memory::MemorySpace memspace) const
-    {
-      using memory::DEVICE;
-      using memory::HOST;
-
-      switch (memspace)
-      {
-      case HOST:
-        if (cpu_updated_[0] == false)
-        {
-          out::error() << "Vector::getData - host data is stale. Perhaps you need to call syncData?\n";
-          return nullptr;
-        }
-        return h_data_;
-      case DEVICE:
-        if (gpu_updated_[0] == false)
-        {
-          out::error() << "Vector::getData - device data is stale. Perhaps you need to call syncData?\n";
-          return nullptr;
-        }
-        return d_data_;
-      default:
-        return nullptr;
-      }
-    }
-
-    /**
-     * @brief Get a pointer to HOST or DEVICE data of a vector in a multivector.
-     *
-     * @param[in] j         - Index of a vector in multivector
-     * @param[in] memspace  - Memory space of the pointer (HOST or DEVICE)
-     *
-     * @return Pointer to the _j_th vector data (HOST or DEVICE).
-     *
-     * @pre `j` < `k_`, i.e., `j` is smaller than the number of vectors.
-     *
-     * @note This function gives you access to the pointer, not to a copy.
-     * If you change the values using the pointer, the vector values will
-     * change too. Call setDataUpdated() to update the staleness flags.
-     */
-    template <typename ScalarT, typename IdxT>
-    ScalarT* Vector<ScalarT, IdxT>::getData(IdxT j, memory::MemorySpace memspace)
-    {
-      using memory::DEVICE;
-      using memory::HOST;
-
-      if (k_ <= j)
-      {
-        out::error() << "Vector::getData - vector index " << j << " out of range, multivector has only " << k_ << " vectors\n";
-        return nullptr;
-      }
-
-      switch (memspace)
-      {
-      case HOST:
-        if (cpu_updated_[j] == false)
-        {
-          out::error() << "Vector::getData - host data for vector " << j << " is stale. Perhaps you need to call syncData?\n";
-          return nullptr;
-        }
-        return &h_data_[j * n_size_];
-      case DEVICE:
-        if (gpu_updated_[j] == false)
-        {
-          out::error() << "Vector::getData - device data for vector " << j << " is stale. Perhaps you need to call syncData?\n";
-          return nullptr;
-        }
-        return &d_data_[j * n_size_];
-      default:
-        return nullptr;
-      }
-    }
-
-    /**
-     * @brief Get a const pointer to HOST or DEVICE data of a vector in a multivector.
-     *
-     * @param[in] j         - Index of a vector in multivector
-     * @param[in] memspace  - Memory space of the pointer (HOST or DEVICE)
-     *
-     * @return Const pointer to the _j_th vector data (HOST or DEVICE).
-     *
-     * @pre `j` < `k_`, i.e., `j` is smaller than the number of vectors.
-     */
-    template <typename ScalarT, typename IdxT>
-    const ScalarT* Vector<ScalarT, IdxT>::getData(IdxT j, memory::MemorySpace memspace) const
-    {
-      using memory::DEVICE;
-      using memory::HOST;
-
-      if (k_ <= j)
-      {
-        out::error() << "Vector::getData - vector index " << j << " out of range, multivector has only " << k_ << " vectors\n";
-        return nullptr;
-      }
-
-      switch (memspace)
-      {
-      case HOST:
-        if (cpu_updated_[j] == false)
-        {
-          out::error() << "Vector::getData - host data for vector " << j << " is stale. Perhaps you need to call syncData?\n";
-          return nullptr;
-        }
-        return &h_data_[j * n_size_];
-      case DEVICE:
-        if (gpu_updated_[j] == false)
-        {
-          out::error() << "Vector::getData - device data for vector " << j << " is stale. Perhaps you need to call syncData?\n";
-          return nullptr;
-        }
-        return &d_data_[j * n_size_];
-      default:
-        return nullptr;
-      }
-    }
-
     /**
      * @brief Sync out of date memory space with the updated one.
      *
@@ -574,19 +378,19 @@ namespace GridKit
     {
       using namespace memory;
 
-      bool all_cpu_updated = cpu_updated_[0];
-      bool all_gpu_updated = gpu_updated_[0];
+      bool all_cpu_updated = getHostUpdated(0);
+      bool all_gpu_updated = getDeviceUpdated(0);
 
       // Verify that all vectors in multivector have the same update status.
       for (IdxT i = 1; i < k_; ++i)
       {
-        if (gpu_updated_[i] != all_gpu_updated)
+        if (getDeviceUpdated(i) != all_gpu_updated)
         {
           out::error() << "Vector::syncData - inconsistent update state across device columns.\n"
                        << "Use syncData(j, memspace) for individual vectors\n";
           return 1;
         }
-        if (cpu_updated_[i] != all_cpu_updated)
+        if (getHostUpdated(i) != all_cpu_updated)
         {
           out::error() << "Vector::syncData - inconsistent update state across host columns.\n"
                        << "Use syncData(j, memspace) for individual vectors\n";
@@ -597,12 +401,12 @@ namespace GridKit
       switch (memspaceDst)
       {
       case DEVICE: // cpu -> gpu
-        if (gpu_updated_[0])
+        if (all_gpu_updated)
         {
           out::error() << "Vector::syncData - device already up to date\n";
           return 1;
         }
-        if (!cpu_updated_[0])
+        if (!all_cpu_updated)
         {
           out::error() << "Vector::syncData - host data is stale, cannot sync to device\n";
           return 1;
@@ -616,12 +420,12 @@ namespace GridKit
         setDeviceUpdated(true);
         break;
       case HOST: // gpu -> cpu
-        if (cpu_updated_[0])
+        if (all_cpu_updated)
         {
           out::error() << "Vector::syncData - host already up to date\n";
           return 1;
         }
-        if (!gpu_updated_[0])
+        if (!all_gpu_updated)
         {
           out::error() << "Vector::syncData - device data is stale, cannot sync to host\n";
           return 1;
@@ -671,12 +475,12 @@ namespace GridKit
       switch (memspaceDst)
       {
       case DEVICE: // cpu->gpu
-        if (gpu_updated_[j])
+        if (getDeviceUpdated(j))
         {
           out::error() << "Vector::syncData - device already up to date\n";
           return 1;
         }
-        if (!cpu_updated_[j])
+        if (!getHostUpdated(j))
         {
           out::error() << "Vector::syncData - host data is stale, cannot sync to device\n";
           return 1;
@@ -687,15 +491,15 @@ namespace GridKit
           return 1;
         }
         mem_.copyArrayHostToDevice(&d_data_[j * n_size_], &h_data_[j * n_size_], n_size_);
-        gpu_updated_[j] = true;
+        setDeviceUpdated(j, true);
         break;
       case HOST: // gpu -> cpu
-        if (cpu_updated_[j])
+        if (getHostUpdated(j))
         {
           out::error() << "Vector::syncData - host already up to date\n";
           return 1;
         }
-        if (!gpu_updated_[j])
+        if (!getDeviceUpdated(j))
         {
           out::error() << "Vector::syncData - device data is stale, cannot sync to host\n";
           return 1;
@@ -706,7 +510,7 @@ namespace GridKit
           return 1;
         }
         mem_.copyArrayDeviceToHost(&h_data_[j * n_size_], &d_data_[j * n_size_], n_size_);
-        cpu_updated_[j] = true;
+        setHostUpdated(j, true);
         break;
       default:
         return 1;
@@ -836,8 +640,8 @@ namespace GridKit
           return 1;
         }
         mem_.setZeroArrayOnHost(&h_data_[j * n_size_], n_size_);
-        cpu_updated_[j] = true;
-        gpu_updated_[j] = false;
+        setHostUpdated(j, true);
+        setDeviceUpdated(j, false);
         break;
       case DEVICE:
         if (d_data_ == nullptr)
@@ -847,8 +651,8 @@ namespace GridKit
         }
         // TODO: We should not need to access raw data in this class
         mem_.setZeroArrayOnDevice(&d_data_[j * n_size_], n_size_);
-        cpu_updated_[j] = false;
-        gpu_updated_[j] = true;
+        setHostUpdated(j, false);
+        setDeviceUpdated(j, true);
         break;
       }
       return 0;
@@ -924,8 +728,8 @@ namespace GridKit
           return 1;
         }
         mem_.setArrayToConstOnHost(&h_data_[n_size_ * j], C, n_size_);
-        cpu_updated_[j] = true;
-        gpu_updated_[j] = false;
+        setHostUpdated(j, true);
+        setDeviceUpdated(j, false);
         break;
       case DEVICE:
         if (d_data_ == nullptr)
@@ -934,8 +738,8 @@ namespace GridKit
           return 1;
         }
         mem_.setArrayToConstOnDevice(&d_data_[n_size_ * j], C, n_size_);
-        cpu_updated_[j] = false;
-        gpu_updated_[j] = true;
+        setHostUpdated(j, false);
+        setDeviceUpdated(j, true);
         break;
       }
       return 0;
@@ -1115,7 +919,7 @@ namespace GridKit
       switch (memspaceSrc)
       {
       case HOST:
-        if (!cpu_updated_[0])
+        if (!getHostUpdated(0))
         {
           out::error() << "Vector::copyToExternal - source data is stale\n";
           return 1;
@@ -1131,7 +935,7 @@ namespace GridKit
         }
         break;
       case DEVICE:
-        if (!gpu_updated_[0])
+        if (!getDeviceUpdated(0))
         {
           out::error() << "Vector::copyToExternal - source data is stale\n";
           return 1;
@@ -1150,22 +954,6 @@ namespace GridKit
       return 0;
     }
 
-    //
-    // Private methods
-    //
-
-    template <typename ScalarT, typename IdxT>
-    void Vector<ScalarT, IdxT>::setHostUpdated(bool is_updated)
-    {
-      std::fill(cpu_updated_, cpu_updated_ + k_, is_updated);
-    }
-
-    template <typename ScalarT, typename IdxT>
-    void Vector<ScalarT, IdxT>::setDeviceUpdated(bool is_updated)
-    {
-      std::fill(gpu_updated_, gpu_updated_ + k_, is_updated);
-    }
-
     template class Vector<double, long int>;
     template class Vector<double, size_t>;
     template class Vector<double, int>;
diff --git a/GridKit/LinearAlgebra/Vector/Vector.hpp b/GridKit/LinearAlgebra/Vector/Vector.hpp
index 823b9149f..bc3a2ff06 100644
--- a/GridKit/LinearAlgebra/Vector/Vector.hpp
+++ b/GridKit/LinearAlgebra/Vector/Vector.hpp
@@ -1,4 +1,6 @@
 #pragma once
+
+#include <cassert>
 #include <string>
 
 #include <GridKit/MemoryUtilities/MemoryUtils.hpp>
@@ -7,6 +9,38 @@ namespace GridKit
 {
   namespace LinearAlgebra
   {
+    namespace detail
+    {
+      // NOTE: this namespace exists to reduce template instantiation + decouple
+      //       code where it isn't necessary
+
+      /**
+       * @brief log failure for host data.
+       */
+      [[gnu::cold, gnu::noinline]]
+      void logHostUnsyncFailure();
+
+      /**
+       * @brief log failure for device data.
+       */
+      [[gnu::cold, gnu::noinline]]
+      void logDeviceUnsyncFailure();
+
+      /**
+       * @brief log failure for bounds check.
+       */
+      template <typename IdxT>
+      [[gnu::cold, gnu::noinline]]
+      void logBoundsCheckFailure(IdxT j, IdxT k);
+
+      /**
+       * @brief log failure for @ref Vector::setDataUpdated bounds check.
+       */
+      template <typename IdxT>
+      [[gnu::cold, gnu::noinline]]
+      void logUpdatedBoundsCheckFailure(IdxT j, IdxT k);
+    } // namespace detail
+
     /**
      * @brief This class implements vectors (dense arrays) and multivectors and
      * some basic utilities (get size, allocate, set data, get data, etc).
@@ -52,17 +86,247 @@ namespace GridKit
                            memory::MemorySpace memspaceIn  = memory::HOST,
                            memory::MemorySpace memspaceOut = memory::HOST);
 
-      ScalarT*       getData(memory::MemorySpace memspace = memory::HOST);
-      ScalarT*       getData(IdxT i, memory::MemorySpace memspace = memory::HOST);
-      const ScalarT* getData(memory::MemorySpace memspace = memory::HOST) const;
-      const ScalarT* getData(IdxT i, memory::MemorySpace memspace = memory::HOST) const;
+      // TODO: in the future, we should consider making the `getData` family of
+      //       methods accept the memory space through a template parameter
+      //       instead to guarantee the optimization implemented here through
+      //       constant propagation
+
+      /**
+       * @brief get a pointer to HOST or DEVICE vector data.
+       *
+       * @param[in] memspace  - Memory space of the pointer (HOST or DEVICE)
+       *
+       * @return pointer to the vector data (HOST or DEVICE). In case of multivectors,
+       * vectors are stored column-wise.
+       *
+       * @note This function gives you access to the pointer, not to a copy.
+       * If you change the values using the pointer, the vector values will
+       * change too. Make sure to use setDataUpdated function to set the update
+       * flags correctly after changing the values.
+       */
+      [[gnu::always_inline]]
+      inline ScalarT* getData(memory::MemorySpace memspace = memory::HOST)
+      {
+        using memory::DEVICE;
+        using memory::HOST;
+
+        switch (memspace)
+        {
+        case HOST:
+          if (!getHostUpdated(0)) [[unlikely]]
+          {
+            detail::logHostUnsyncFailure();
+            return nullptr;
+          }
+          return h_data_;
+        case DEVICE:
+          if (!getDeviceUpdated(0)) [[unlikely]]
+          {
+            detail::logDeviceUnsyncFailure();
+            return nullptr;
+          }
+          return d_data_;
+        default:
+          return nullptr;
+        }
+      }
+
+      /**
+       * @brief Get a pointer to HOST or DEVICE data of a vector in a multivector.
+       *
+       * @param[in] j         - Index of a vector in multivector
+       * @param[in] memspace  - Memory space of the pointer (HOST or DEVICE)
+       *
+       * @return Pointer to the _j_th vector data (HOST or DEVICE).
+       *
+       * @pre `j` < `k_`, i.e., `j` is smaller than the number of vectors.
+       *
+       * @note This function gives you access to the pointer, not to a copy.
+       * If you change the values using the pointer, the vector values will
+       * change too. Call setDataUpdated() to update the staleness flags.
+       */
+      [[gnu::always_inline]]
+      ScalarT* getData(IdxT j, memory::MemorySpace memspace = memory::HOST)
+      {
+        using memory::DEVICE;
+        using memory::HOST;
+
+        if (k_ <= j) [[unlikely]]
+        {
+          detail::logBoundsCheckFailure<IdxT>(j, k_);
+          return nullptr;
+        }
+
+        switch (memspace)
+        {
+        case HOST:
+          if (!getHostUpdated(j)) [[unlikely]]
+          {
+            detail::logHostUnsyncFailure();
+            return nullptr;
+          }
+          return &h_data_[j * n_size_];
+        case DEVICE:
+          if (!getDeviceUpdated(j)) [[unlikely]]
+          {
+            detail::logDeviceUnsyncFailure();
+            return nullptr;
+          }
+          return &d_data_[j * n_size_];
+        default:
+          return nullptr;
+        }
+      }
+
+      /**
+       * @brief get a pointer to HOST or DEVICE vector data.
+       *
+       * @param[in] memspace  - Memory space of the pointer (HOST or DEVICE)
+       *
+       * @return pointer to the vector data (HOST or DEVICE). In case of multivectors,
+       * vectors are stored column-wise.
+       */
+      [[gnu::always_inline]]
+      const ScalarT* getData(memory::MemorySpace memspace = memory::HOST) const
+      {
+        using memory::DEVICE;
+        using memory::HOST;
+
+        switch (memspace)
+        {
+        case HOST:
+          if (!getHostUpdated(0)) [[unlikely]]
+          {
+            detail::logHostUnsyncFailure();
+            return nullptr;
+          }
+          return h_data_;
+        case DEVICE:
+          if (!getDeviceUpdated(0)) [[unlikely]]
+          {
+            detail::logDeviceUnsyncFailure();
+            return nullptr;
+          }
+          return d_data_;
+        default:
+          return nullptr;
+        }
+      }
+
+      /**
+       * @brief Get a const pointer to HOST or DEVICE data of a vector in a multivector.
+       *
+       * @param[in] j         - Index of a vector in multivector
+       * @param[in] memspace  - Memory space of the pointer (HOST or DEVICE)
+       *
+       * @return Const pointer to the _j_th vector data (HOST or DEVICE).
+       *
+       * @pre `j` < `k_`, i.e., `j` is smaller than the number of vectors.
+       */
+      [[gnu::always_inline]]
+      const ScalarT* getData(IdxT j, memory::MemorySpace memspace = memory::HOST) const
+      {
+        using memory::DEVICE;
+        using memory::HOST;
+
+        if (k_ <= j) [[unlikely]]
+        {
+          detail::logBoundsCheckFailure<IdxT>(j, k_);
+          return nullptr;
+        }
+
+        switch (memspace)
+        {
+        case HOST:
+          if (!getHostUpdated(j)) [[unlikely]]
+          {
+            detail::logHostUnsyncFailure();
+            return nullptr;
+          }
+          return &h_data_[j * n_size_];
+        case DEVICE:
+          if (!getDeviceUpdated(j)) [[unlikely]]
+          {
+            detail::logDeviceUnsyncFailure();
+            return nullptr;
+          }
+          return &d_data_[j * n_size_];
+        default:
+          return nullptr;
+        }
+      }
 
       IdxT getCapacity() const;
       IdxT getSize() const;
       IdxT getNumVectors() const;
 
-      int setDataUpdated(memory::MemorySpace memspace = memory::HOST);
-      int setDataUpdated(IdxT j, memory::MemorySpace memspace = memory::HOST);
+      /**
+       * @brief Set the flag to indicate that the data (HOST or DEVICE) has been
+       * updated.
+       *
+       * Use this function if you update vector elements by accessing the raw data
+       * pointer.
+       *
+       * @param[in] memspace - Memory space (HOST or DEVICE)
+       *
+       * @warning This is an expert level method. Use only if you know what
+       * you are doing.
+       */
+      [[gnu::always_inline]]
+      int setDataUpdated(memory::MemorySpace memspace = memory::HOST)
+      {
+        using namespace memory;
+        switch (memspace)
+        {
+        case HOST:
+          setHostUpdated(true);
+          setDeviceUpdated(false);
+          break;
+        case DEVICE:
+          setHostUpdated(false);
+          setDeviceUpdated(true);
+          break;
+        }
+        return 0;
+      }
+
+      /**
+       * @brief Set the flag to indicate that the data (HOST or DEVICE) for
+       * vector `j` in the multivector has been updated.
+       *
+       * Use this function if you update vector elements by accessing the raw data
+       * pointer.
+       *
+       * @param[in] memspace - Memory space (HOST or DEVICE)
+       *
+       * @warning This is an expert level method. Use only if you know what
+       * you are doing.
+       */
+      [[gnu::always_inline]]
+      int setDataUpdated(IdxT j, memory::MemorySpace memspace = memory::HOST)
+      {
+        using namespace memory;
+
+        if (k_ <= j) [[unlikely]]
+        {
+          detail::logUpdatedBoundsCheckFailure(j, k_);
+          return 1;
+        }
+
+        switch (memspace)
+        {
+        case HOST:
+          setHostUpdated(j, true);
+          setDeviceUpdated(j, false);
+          break;
+        case DEVICE:
+          setDeviceUpdated(j, true);
+          setHostUpdated(j, false);
+          break;
+        }
+        return 0;
+      }
+
       int setData(ScalarT* data, memory::MemorySpace memspace = memory::HOST);
       int setData(ScalarT* data, IdxT size, memory::MemorySpace memspace = memory::HOST);
       int allocate(memory::MemorySpace memspace = memory::HOST);
@@ -82,18 +346,107 @@ namespace GridKit
                          memory::MemorySpace memspaceDst = memory::HOST);
 
     private:
-      void setHostUpdated(bool is_updated);
-      void setDeviceUpdated(bool is_updated);
-
-      IdxT     n_capacity_{0};        ///< vector capacity
-      IdxT     k_{0};                 ///< number of vectors in multivector
-      IdxT     n_size_{0};            ///< actual size of the vector
-      ScalarT* d_data_{nullptr};      ///< DEVICE data array
-      ScalarT* h_data_{nullptr};      ///< HOST data array
-      bool*    gpu_updated_{nullptr}; ///< DEVICE data flags (updated or not)
-      bool*    cpu_updated_{nullptr}; ///< HOST data flags (updated or not)
-
-      bool owns_gpu_data_{true}; ///< data owneship flag for DEVICE data
+      [[gnu::always_inline]]
+      void setHostUpdated(bool is_updated)
+      {
+        if (k_ <= 1) [[likely]]
+        {
+          cpu_updated_spec_ = is_updated;
+          return;
+        }
+
+        std::fill(cpu_updated_, cpu_updated_ + k_, is_updated);
+      }
+
+      [[gnu::always_inline]]
+      void setHostUpdated(IdxT j, bool is_updated)
+      {
+        if (k_ <= 1)
+        {
+          assert(j == 0);
+          cpu_updated_spec_ = is_updated;
+          return;
+        }
+
+        cpu_updated_[j] = is_updated;
+      }
+
+      [[gnu::always_inline]]
+      void setDeviceUpdated(bool is_updated)
+      {
+        if (k_ <= 1) [[likely]]
+        {
+          gpu_updated_spec_ = is_updated;
+          return;
+        }
+
+        std::fill(gpu_updated_, gpu_updated_ + k_, is_updated);
+      }
+
+      [[gnu::always_inline]]
+      void setDeviceUpdated(IdxT j, bool is_updated)
+      {
+        if (k_ <= 1)
+        {
+          assert(j == 0);
+          gpu_updated_spec_ = is_updated;
+          return;
+        }
+
+        gpu_updated_[j] = is_updated;
+      }
+
+      [[gnu::always_inline]]
+      bool getHostUpdated(IdxT j) const
+      {
+        if (k_ <= 1) [[likely]]
+        {
+          assert(j == 0);
+          return cpu_updated_spec_;
+        }
+
+        return cpu_updated_[j];
+      }
+
+      [[gnu::always_inline]]
+      bool getDeviceUpdated(IdxT j) const
+      {
+        if (k_ <= 1) [[likely]]
+        {
+          assert(j == 0);
+          return gpu_updated_spec_;
+        }
+
+        return gpu_updated_[j];
+      }
+
+      IdxT     n_capacity_{0};   ///< vector capacity
+      IdxT     k_{0};            ///< number of vectors in multivector
+      IdxT     n_size_{0};       ///< actual size of the vector
+      ScalarT* d_data_{nullptr}; ///< DEVICE data array
+      ScalarT* h_data_{nullptr}; ///< HOST data array
+
+      /// device data flags
+      union
+      {
+        /// when k > 1, store the flags on the heap
+        bool* gpu_updated_{nullptr};
+
+        /// when k == 1, store the flag in-place for locality
+        bool gpu_updated_spec_;
+      };
+
+      /// host data flags
+      union
+      {
+        /// when k > 1, store the flags on the heap
+        bool* cpu_updated_{nullptr};
+
+        /// when k == 1, store the flag in-place for locality
+        bool cpu_updated_spec_;
+      };
+
+      bool owns_gpu_data_{true}; ///< data ownership flag for DEVICE data
       bool owns_cpu_data_{true}; ///< data ownership flag for HOST data
 
       MemoryManager mem_; ///< Device memory manager object
diff --git a/GridKit/Model/PhasorDynamics/Branch/BranchDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Branch/BranchDependencyTracking.cpp
index 547e8526c..ad068b8ab 100644
--- a/GridKit/Model/PhasorDynamics/Branch/BranchDependencyTracking.cpp
+++ b/GridKit/Model/PhasorDynamics/Branch/BranchDependencyTracking.cpp
@@ -22,9 +22,6 @@ namespace GridKit
     template <typename scalar_type, typename index_type>
     int Branch<scalar_type, index_type>::evaluateJacobian()
     {
-      Log::misc() << "Evaluate DependencyTracking Jacobian for Branch...\n";
-      Log::misc() << "Jacobian evaluation is experimental!\n";
-
       return 0;
     }
 
diff --git a/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp b/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp
index 9094a60ac..c742968b6 100644
--- a/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp
@@ -20,9 +20,6 @@ namespace GridKit
     template <typename scalar_type, typename index_type>
     int Branch<scalar_type, index_type>::evaluateJacobian()
     {
-      Log::misc() << "Evaluate Jacobian for Branch..." << std::endl;
-      Log::misc() << "Jacobian evaluation is experimental!" << std::endl;
-
       if (J_rows_buffer_ == nullptr)
       {
         // Reserve space for the dense blocks.
diff --git a/GridKit/Model/PhasorDynamics/Bus/BusDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Bus/BusDependencyTracking.cpp
index a24ca20b8..c8ae5d8c4 100644
--- a/GridKit/Model/PhasorDynamics/Bus/BusDependencyTracking.cpp
+++ b/GridKit/Model/PhasorDynamics/Bus/BusDependencyTracking.cpp
@@ -24,9 +24,6 @@ namespace GridKit
     template <typename scalar_type, typename index_type>
     int Bus<scalar_type, index_type>::evaluateJacobian()
     {
-      Log::misc() << "Evaluate DependencyTracking Jacobian for Bus...\n";
-      Log::misc() << "Jacobian evaluation is not implemented!\n";
-
       return 0;
     }
 
diff --git a/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp b/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp
index 4af2d5a35..25230d60c 100644
--- a/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp
@@ -22,9 +22,6 @@ namespace GridKit
     template <typename scalar_type, typename index_type>
     int Bus<scalar_type, index_type>::evaluateJacobian()
     {
-      Log::misc() << "Evaluate Jacobian for Bus..." << std::endl;
-      Log::misc() << "Jacobian evaluation is experimental!" << std::endl;
-
       if (J_rows_buffer_ == nullptr)
       {
         J_rows_buffer_ = new IdxT[4];
diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFaultDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/BusFault/BusFaultDependencyTracking.cpp
index b50cc68d9..c9d27f538 100644
--- a/GridKit/Model/PhasorDynamics/BusFault/BusFaultDependencyTracking.cpp
+++ b/GridKit/Model/PhasorDynamics/BusFault/BusFaultDependencyTracking.cpp
@@ -23,9 +23,6 @@ namespace GridKit
     template <class scalar_type, typename index_type>
     int BusFault<scalar_type, index_type>::evaluateJacobian()
     {
-      Log::misc() << "Evaluate DependencyTracking Jacobian for BusFault...\n";
-      Log::misc() << "Jacobian evaluation is experimental!\n";
-
       this->constructCsr();
 
       return 0;
diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp b/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp
index 359e79f23..ba88ca209 100644
--- a/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp
@@ -20,9 +20,6 @@ namespace GridKit
     template <class scalar_type, typename index_type>
     int BusFault<scalar_type, index_type>::evaluateJacobian()
     {
-      Log::misc() << "Evaluate Jacobian for BusFault..." << std::endl;
-      Log::misc() << "Jacobian evaluation is experimental!" << std::endl;
-
       if (J_rows_buffer_ == nullptr)
       {
         // Reserve space for the dense blocks.
diff --git a/GridKit/Model/PhasorDynamics/Component.hpp b/GridKit/Model/PhasorDynamics/Component.hpp
index 0c1529923..a4bf74c7b 100644
--- a/GridKit/Model/PhasorDynamics/Component.hpp
+++ b/GridKit/Model/PhasorDynamics/Component.hpp
@@ -264,7 +264,7 @@ namespace GridKit
       /**
        * @brief CSR construction dispatch depending on ScalarT
        *
-       * @note Currently only used for testing, particularly for comparing 
+       * @note Currently only used for testing, particularly for comparing
        *       Enzyme and DependencyTracking Jacobians.
        */
       int constructCsr()
@@ -340,7 +340,7 @@ namespace GridKit
       /**
        * @brief CSR construction from COO.
        *
-       * @note Currently only used for testing, particularly for comparing 
+       * @note Currently only used for testing, particularly for comparing
        *       Enzyme and DependencyTracking Jacobians.
        *
        * @todo The matrix is only computed on the first call, and the data is stale on subsequent calls.
@@ -374,17 +374,17 @@ namespace GridKit
       /**
        * @brief CSR construction from Dependency maps.
        *
-       * This merges the y and yp dependencies back to the same location in the Jacobian. 
+       * This merges the y and yp dependencies back to the same location in the Jacobian.
        * See \ref initializeDependencyTrackingVariableNumbers() for the initial even/odd split.
-       * The DependencyMap manipulations are expected to be a bottleneck, so this is not intended 
-       * for simulations. 
+       * The DependencyMap manipulations are expected to be a bottleneck, so this is not intended
+       * for simulations.
        *
-       * @note Currently only used for testing, particularly for comparing 
-       *       Enzyme and DependencyTracking Jacobians. Will be used to benchmark the 
+       * @note Currently only used for testing, particularly for comparing
+       *       Enzyme and DependencyTracking Jacobians. Will be used to benchmark the
        *       performance of the two methods in the near future.
        *
-       * @note For dependency maps, there is no functional difference in Jacobian construction 
-       *       between components and systems. However, the resulting component-level CSR will 
+       * @note For dependency maps, there is no functional difference in Jacobian construction
+       *       between components and systems. However, the resulting component-level CSR will
        *       be invalid if tracked indices are greater than the local `size_`, and should not
        *       be used for simulation. Systems do not rely on the component-level results.
        */
diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbDependencyTracking.cpp
index b7dd60c5f..b1841f428 100644
--- a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbDependencyTracking.cpp
+++ b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbDependencyTracking.cpp
@@ -23,9 +23,6 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int Reecb<scalar_type, index_type>::evaluateJacobian()
       {
-        Log::misc() << "Evaluate DependencyTracking Jacobian for Reecb...\n";
-        Log::misc() << "Jacobian evaluation is experimental!\n";
-
         this->constructCsr();
 
         return 0;
diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbEnzyme.cpp b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbEnzyme.cpp
index 682b18700..f474aa176 100644
--- a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbEnzyme.cpp
@@ -30,9 +30,6 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int Reecb<scalar_type, index_type>::evaluateJacobian()
       {
-        Log::misc() << "Evaluate Jacobian for Reecb...\n";
-        Log::misc() << "Jacobian evaluation is experimental!\n";
-
         if (J_rows_buffer_ == nullptr)
         {
           const auto size        = static_cast<size_t>(size_);
diff --git a/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaDependencyTracking.cpp
index 185097dda..8b464d191 100644
--- a/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaDependencyTracking.cpp
+++ b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaDependencyTracking.cpp
@@ -23,9 +23,6 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int Repca<scalar_type, index_type>::evaluateJacobian()
       {
-        Log::misc() << "Evaluate DependencyTracking Jacobian for Repca...\n";
-        Log::misc() << "Jacobian evaluation is experimental!\n";
-
         this->constructCsr();
 
         return 0;
diff --git a/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaEnzyme.cpp b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaEnzyme.cpp
index 1b0242296..b6ee25958 100644
--- a/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaEnzyme.cpp
@@ -30,9 +30,6 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int Repca<scalar_type, index_type>::evaluateJacobian()
       {
-        Log::misc() << "Evaluate Jacobian for Repca...\n";
-        Log::misc() << "Jacobian evaluation is experimental!\n";
-
         if (J_rows_buffer_ == nullptr)
         {
           const auto size        = static_cast<size_t>(size_);
diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp
index 6945845d9..39f258514 100644
--- a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp
+++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp
@@ -23,9 +23,6 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int Regca<scalar_type, index_type>::evaluateJacobian()
       {
-        Log::misc() << "Evaluate DependencyTracking Jacobian for Regca...\n";
-        Log::misc() << "Jacobian evaluation is experimental!\n";
-
         this->constructCsr();
 
         return 0;
diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp
index 935c00320..c68327fcd 100644
--- a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp
@@ -22,9 +22,6 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int Regca<scalar_type, index_type>::evaluateJacobian()
       {
-        Log::misc() << "Evaluate Jacobian for Regca..." << std::endl;
-        Log::misc() << "Jacobian evaluation is experimental!" << std::endl;
-
         if (J_rows_buffer_ == nullptr)
         {
           // Reserve space for the dense blocks. Enzyme keeps only structural
diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aDependencyTracking.cpp
index 7fd338323..02470c673 100644
--- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aDependencyTracking.cpp
+++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aDependencyTracking.cpp
@@ -23,9 +23,6 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int Esdc1a<scalar_type, index_type>::evaluateJacobian()
       {
-        Log::misc() << "Evaluate DependencyTracking Jacobian for Esdc1a...\n";
-        Log::misc() << "Jacobian evaluation is experimental!\n";
-
         this->constructCsr();
 
         return 0;
diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aEnzyme.cpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aEnzyme.cpp
index d3ea94be6..45c38313b 100644
--- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aEnzyme.cpp
@@ -29,9 +29,6 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int Esdc1a<scalar_type, index_type>::evaluateJacobian()
       {
-        Log::misc() << "Evaluate Jacobian for Esdc1a..." << std::endl;
-        Log::misc() << "Jacobian evaluation is experimental!" << std::endl;
-
         if (J_rows_buffer_ == nullptr)
         {
           auto size        = static_cast<size_t>(size_);
diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1DependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1DependencyTracking.cpp
index 43df43a0d..b0a8e9457 100644
--- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1DependencyTracking.cpp
+++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1DependencyTracking.cpp
@@ -28,9 +28,6 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int Ieeet1<scalar_type, index_type>::evaluateJacobian()
       {
-        Log::misc() << "Evaluate DependencyTracking Jacobian for Ieeet1...\n";
-        Log::misc() << "Jacobian evaluation is experimental!\n";
-
         this->constructCsr();
 
         return 0;
diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp
index adf3948a0..1ece23a46 100644
--- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp
@@ -22,9 +22,6 @@ namespace GridKit
       template <class scalar_type, typename index_type>
       int Ieeet1<scalar_type, index_type>::evaluateJacobian()
       {
-        Log::misc() << "Evaluate Jacobian for Ieeet1..." << std::endl;
-        Log::misc() << "Jacobian evaluation is experimental!" << std::endl;
-
         if (J_rows_buffer_ == nullptr)
         {
           // Reserve space for the dense blocks.
diff --git a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiDependencyTracking.cpp
index e79ad9715..1e01831ef 100644
--- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiDependencyTracking.cpp
+++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiDependencyTracking.cpp
@@ -25,9 +25,6 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int SexsPti<scalar_type, index_type>::evaluateJacobian()
       {
-        Log::misc() << "Evaluate DependencyTracking Jacobian for SexsPti...\n";
-        Log::misc() << "Jacobian evaluation is experimental!\n";
-
         this->constructCsr();
 
         return 0;
diff --git a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiEnzyme.cpp b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiEnzyme.cpp
index a15c31c10..2687ef804 100644
--- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiEnzyme.cpp
@@ -17,9 +17,6 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int SexsPti<scalar_type, index_type>::evaluateJacobian()
       {
-        Log::misc() << "Evaluate Jacobian for SexsPti..." << std::endl;
-        Log::misc() << "Jacobian evaluation is experimental!" << std::endl;
-
         if (J_rows_buffer_ == nullptr)
         {
           // Reserve space for the dense blocks.
diff --git a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiDependencyTracking.cpp
index 6dde71a2f..4d95bfbce 100644
--- a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiDependencyTracking.cpp
+++ b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiDependencyTracking.cpp
@@ -25,9 +25,6 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int GastPti<scalar_type, index_type>::evaluateJacobian()
       {
-        Log::misc() << "Evaluate DependencyTracking Jacobian for GastPti...\n";
-        Log::misc() << "Jacobian evaluation is experimental!\n";
-
         this->constructCsr();
 
         return 0;
diff --git a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiEnzyme.cpp b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiEnzyme.cpp
index c41d5a7ff..829812a75 100644
--- a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiEnzyme.cpp
@@ -29,9 +29,6 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int GastPti<scalar_type, index_type>::evaluateJacobian()
       {
-        Log::misc() << "Evaluate Jacobian for GastPti..." << std::endl;
-        Log::misc() << "Jacobian evaluation is experimental!" << std::endl;
-
         if (J_rows_buffer_ == nullptr)
         {
           auto size        = static_cast<size_t>(size_);
diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovDependencyTracking.cpp
index 406faf306..adb4d58fd 100644
--- a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovDependencyTracking.cpp
+++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovDependencyTracking.cpp
@@ -25,9 +25,6 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int Hygov<scalar_type, index_type>::evaluateJacobian()
       {
-        Log::misc() << "Evaluate DependencyTracking Jacobian for Hygov...\n";
-        Log::misc() << "Jacobian evaluation is experimental!\n";
-
         this->constructCsr();
 
         return 0;
diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp
index af4cb2c90..acd5f3db6 100644
--- a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp
@@ -17,9 +17,6 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int Hygov<scalar_type, index_type>::evaluateJacobian()
       {
-        Log::misc() << "Evaluate Jacobian for Hygov..." << std::endl;
-        Log::misc() << "Jacobian evaluation is experimental!" << std::endl;
-
         if (J_rows_buffer_ == nullptr)
         {
           auto size        = static_cast<size_t>(size_);
diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp
index 39e730e6b..dd78ad007 100644
--- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Enzyme.cpp
@@ -22,9 +22,6 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int Tgov1<scalar_type, index_type>::evaluateJacobian()
       {
-        Log::misc() << "Evaluate Jacobian for Tgov1..." << std::endl;
-        Log::misc() << "Jacobian evaluation is experimental!" << std::endl;
-
         if (J_rows_buffer_ == nullptr)
         {
           // Reserve space for the dense blocks.
diff --git a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZDependencyTracking.cpp
index b02b683a2..033d70586 100644
--- a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZDependencyTracking.cpp
+++ b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZDependencyTracking.cpp
@@ -18,9 +18,6 @@ namespace GridKit
     template <typename scalar_type, typename index_type>
     int LoadZ<scalar_type, index_type>::evaluateJacobian()
     {
-      Log::misc() << "Evaluate DependencyTracking Jacobian for LoadZ...\n";
-      Log::misc() << "Jacobian evaluation is experimental!\n";
-
       this->constructCsr();
 
       return 0;
diff --git a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZEnzyme.cpp b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZEnzyme.cpp
index 84a03966d..965390e22 100644
--- a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZEnzyme.cpp
@@ -20,9 +20,6 @@ namespace GridKit
     template <typename scalar_type, typename index_type>
     int LoadZ<scalar_type, index_type>::evaluateJacobian()
     {
-      Log::misc() << "Evaluate Jacobian for LoadZ..." << std::endl;
-      Log::misc() << "Jacobian evaluation is experimental!" << std::endl;
-
       if (J_rows_buffer_ == nullptr)
       {
         // Reserve space for the dense blocks.
diff --git a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPDependencyTracking.cpp
index 2c84ee6b1..66cd4fbba 100644
--- a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPDependencyTracking.cpp
+++ b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPDependencyTracking.cpp
@@ -18,9 +18,6 @@ namespace GridKit
     template <typename scalar_type, typename index_type>
     int LoadZIP<scalar_type, index_type>::evaluateJacobian()
     {
-      Log::misc() << "Evaluate DependencyTracking Jacobian for LoadZIP...\n";
-      Log::misc() << "Jacobian evaluation is experimental!\n";
-
       this->constructCsr();
 
       return 0;
diff --git a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPEnzyme.cpp b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPEnzyme.cpp
index 4e25242f7..6465b08a5 100644
--- a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPEnzyme.cpp
@@ -15,9 +15,6 @@ namespace GridKit
     template <typename scalar_type, typename index_type>
     int LoadZIP<scalar_type, index_type>::evaluateJacobian()
     {
-      Log::misc() << "Evaluate Jacobian for LoadZIP..." << std::endl;
-      Log::misc() << "Jacobian evaluation is experimental!" << std::endl;
-
       if (J_rows_buffer_ == nullptr)
       {
         // Reserve space for the dense blocks.
diff --git a/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.cpp b/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.cpp
index d0310c67d..0502a92eb 100644
--- a/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.cpp
+++ b/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.cpp
@@ -1,7 +1,7 @@
 /**
  * @file SignalNode model implementation.
  */
-#include "SignalNodeImpl.hpp"
+#include "SignalNode.hpp"
 
 namespace GridKit
 {
diff --git a/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp b/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp
index fa745fac4..7c9e8118f 100644
--- a/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp
+++ b/GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp
@@ -1,17 +1,11 @@
 #pragma once
 
+#include <cassert>
+
 #include <GridKit/Constants.hpp>
+#include <GridKit/Model/PhasorDynamics/SignalNode/SignalNodeData.hpp>
 #include <GridKit/ScalarTraits.hpp>
 
-namespace GridKit
-{
-  namespace PhasorDynamics
-  {
-    template <typename real_type, typename index_type>
-    struct SignalNodeData;
-  } // namespace PhasorDynamics
-} // namespace GridKit
-
 namespace GridKit
 {
   namespace PhasorDynamics
@@ -39,19 +33,69 @@ namespace GridKit
       using IdxT    = index_type;
       using RealT   = typename GridKit::ScalarTraits<ScalarT>::RealT;
 
-      SignalNode();
-      SignalNode(const SignalNodeData<RealT, IdxT>& data);
+      SignalNode()
+      {
+      }
+
+      SignalNode(const SignalNodeData<RealT, IdxT>& data)
+        : signal_id_(data.signal_id)
+      {
+      }
 
       virtual ~SignalNode() = default;
 
-      IdxT    signalId() const noexcept;
-      void    setAssigned() noexcept;
-      bool    assigned() const noexcept;
-      void    link(ScalarT* signal_in, IdxT* global_index) noexcept;
-      bool    linked() const noexcept;
-      ScalarT read() const noexcept;
-      IdxT    getVariableIndex() const noexcept;
-      void    init(ScalarT signal_in) noexcept;
+      [[gnu::always_inline]]
+      IdxT signalId() const noexcept
+      {
+        return signal_id_;
+      }
+
+      [[gnu::always_inline]]
+      void setAssigned() noexcept
+      {
+        assert(!assigned_);
+        assigned_ = true;
+      }
+
+      [[gnu::always_inline]]
+      bool assigned() const noexcept
+      {
+        return assigned_;
+      }
+
+      [[gnu::always_inline]]
+      void link(ScalarT* signal_in, IdxT* global_index) noexcept
+      {
+        signal_         = signal_in;
+        variable_index_ = global_index;
+      }
+
+      [[gnu::always_inline]]
+      bool linked() const noexcept
+      {
+        return (signal_) && (variable_index_);
+      }
+
+      [[gnu::always_inline]]
+      ScalarT read() const noexcept
+      {
+        assert(signal_);
+        return *signal_;
+      }
+
+      [[gnu::always_inline]]
+      IdxT getVariableIndex() const noexcept
+      {
+        assert(variable_index_);
+        return *variable_index_;
+      }
+
+      [[gnu::always_inline]]
+      void init(ScalarT signal_in) noexcept
+      {
+        assert(signal_);
+        *signal_ = signal_in;
+      }
 
     private:
       ScalarT* signal_{nullptr};
diff --git a/GridKit/Model/PhasorDynamics/SignalNode/SignalNodeDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SignalNode/SignalNodeDependencyTracking.cpp
index 33dc3810d..65401ecc8 100644
--- a/GridKit/Model/PhasorDynamics/SignalNode/SignalNodeDependencyTracking.cpp
+++ b/GridKit/Model/PhasorDynamics/SignalNode/SignalNodeDependencyTracking.cpp
@@ -4,7 +4,7 @@
 
 #include <GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp>
 
-#include "SignalNodeImpl.hpp"
+#include "SignalNode.hpp"
 
 namespace GridKit
 {
diff --git a/GridKit/Model/PhasorDynamics/SignalNode/SignalNodeImpl.hpp b/GridKit/Model/PhasorDynamics/SignalNode/SignalNodeImpl.hpp
deleted file mode 100644
index cd32c84ed..000000000
--- a/GridKit/Model/PhasorDynamics/SignalNode/SignalNodeImpl.hpp
+++ /dev/null
@@ -1,78 +0,0 @@
-/**
- * @file SignalNode model implementation.
- */
-
-#include <cassert>
-
-#include <GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp>
-#include <GridKit/Model/PhasorDynamics/SignalNode/SignalNodeData.hpp>
-
-namespace GridKit
-{
-  namespace PhasorDynamics
-  {
-    template <typename scalar_type, typename index_type>
-    SignalNode<scalar_type, index_type>::SignalNode()
-    {
-    }
-
-    template <typename scalar_type, typename index_type>
-    SignalNode<scalar_type, index_type>::SignalNode(const SignalNodeData<RealT, IdxT>& data)
-      : signal_id_(data.signal_id)
-    {
-    }
-
-    template <typename scalar_type, typename index_type>
-    index_type SignalNode<scalar_type, index_type>::signalId() const noexcept
-    {
-      return signal_id_;
-    }
-
-    template <typename scalar_type, typename index_type>
-    void SignalNode<scalar_type, index_type>::setAssigned() noexcept
-    {
-      assert(!assigned_);
-      assigned_ = true;
-    }
-
-    template <typename scalar_type, typename index_type>
-    bool SignalNode<scalar_type, index_type>::assigned() const noexcept
-    {
-      return assigned_;
-    }
-
-    template <typename scalar_type, typename index_type>
-    void SignalNode<scalar_type, index_type>::link(ScalarT* signal, IdxT* variable_index) noexcept
-    {
-      signal_         = signal;
-      variable_index_ = variable_index;
-    }
-
-    template <typename scalar_type, typename index_type>
-    bool SignalNode<scalar_type, index_type>::linked() const noexcept
-    {
-      return (signal_) && (variable_index_);
-    }
-
-    template <typename scalar_type, typename index_type>
-    scalar_type SignalNode<scalar_type, index_type>::read() const noexcept
-    {
-      assert(signal_);
-      return *signal_;
-    }
-
-    template <typename scalar_type, typename index_type>
-    index_type SignalNode<scalar_type, index_type>::getVariableIndex() const noexcept
-    {
-      assert(variable_index_);
-      return *variable_index_;
-    }
-
-    template <typename scalar_type, typename index_type>
-    void SignalNode<scalar_type, index_type>::init(ScalarT signal) noexcept
-    {
-      assert(signal_);
-      *signal_ = signal;
-    }
-  } // namespace PhasorDynamics
-} // namespace GridKit
diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestDependencyTracking.cpp
index eed3e9bd9..376b1fe21 100644
--- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestDependencyTracking.cpp
+++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestDependencyTracking.cpp
@@ -25,9 +25,6 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int Ieeest<scalar_type, index_type>::evaluateJacobian()
       {
-        Log::misc() << "Evaluate DependencyTracking Jacobian for Ieeest...\n";
-        Log::misc() << "Jacobian evaluation is experimental!\n";
-
         this->constructCsr();
 
         return 0;
diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp
index c814afe3e..eabaed17d 100644
--- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp
@@ -24,9 +24,6 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int Ieeest<scalar_type, index_type>::evaluateJacobian()
       {
-        Log::misc() << "Evaluate Jacobian for Ieeest..." << std::endl;
-        Log::misc() << "Jacobian evaluation is experimental!" << std::endl;
-
         if (J_rows_buffer_ == nullptr)
         {
           // Reserve space for the dense blocks.
diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouDependencyTracking.cpp
index fdf9b1805..2f3bd55be 100644
--- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouDependencyTracking.cpp
+++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouDependencyTracking.cpp
@@ -17,9 +17,6 @@ namespace GridKit
     template <typename scalar_type, typename index_type>
     int Genrou<scalar_type, index_type>::evaluateJacobian()
     {
-      Log::misc() << "Evaluate DependencyTracking Jacobian for Genrou...\n";
-      Log::misc() << "Jacobian evaluation is experimental!\n";
-
       this->constructCsr();
 
       return 0;
diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouEnzyme.cpp
index 242d92eaf..5d473db8f 100644
--- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouEnzyme.cpp
@@ -20,9 +20,6 @@ namespace GridKit
     template <typename scalar_type, typename index_type>
     int Genrou<scalar_type, index_type>::evaluateJacobian()
     {
-      Log::misc() << "Evaluate Jacobian for Genrou..." << std::endl;
-      Log::misc() << "Jacobian evaluation is experimental!" << std::endl;
-
       if (J_rows_buffer_ == nullptr)
       {
         // Reserve space for the dense blocks.
diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalDependencyTracking.cpp
index f366f8029..bc492da06 100644
--- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalDependencyTracking.cpp
+++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalDependencyTracking.cpp
@@ -17,9 +17,6 @@ namespace GridKit
     template <typename scalar_type, typename index_type>
     int Gensal<scalar_type, index_type>::evaluateJacobian()
     {
-      Log::misc() << "Evaluate DependencyTracking Jacobian for Gensal...\n";
-      Log::misc() << "Jacobian evaluation is experimental!\n";
-
       this->constructCsr();
 
       return 0;
diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalEnzyme.cpp
index 810a725c8..354c8428c 100644
--- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalEnzyme.cpp
@@ -20,9 +20,6 @@ namespace GridKit
     template <typename scalar_type, typename index_type>
     int Gensal<scalar_type, index_type>::evaluateJacobian()
     {
-      Log::misc() << "Evaluate Jacobian for Gensal..." << std::endl;
-      Log::misc() << "Jacobian evaluation is experimental!" << std::endl;
-
       if (J_rows_buffer_ == nullptr)
       {
         // Reserve space for the dense blocks.
diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalDependencyTracking.cpp
index e4eebae9e..3b31c9fd0 100644
--- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalDependencyTracking.cpp
+++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalDependencyTracking.cpp
@@ -18,9 +18,6 @@ namespace GridKit
     template <typename scalar_type, typename index_type>
     int GenClassical<scalar_type, index_type>::evaluateJacobian()
     {
-      Log::misc() << "Evaluate DependencyTracking Jacobian for GenClassical...\n";
-      Log::misc() << "Jacobian evaluation is experimental!\n";
-
       this->constructCsr();
 
       return 0;
diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp
index 0f8ea63ec..72475bd90 100644
--- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp
@@ -20,9 +20,6 @@ namespace GridKit
     template <typename scalar_type, typename index_type>
     int GenClassical<scalar_type, index_type>::evaluateJacobian()
     {
-      Log::misc() << "Evaluate Jacobian for GenClassical..." << std::endl;
-      Log::misc() << "Jacobian evaluation is experimental!" << std::endl;
-
       if (J_rows_buffer_ == nullptr)
       {
         // Reserve space for the dense blocks.
diff --git a/GridKit/Solver/Dynamic/Ida.cpp b/GridKit/Solver/Dynamic/Ida.cpp
index 1b3de5658..03a5e4dd0 100644
--- a/GridKit/Solver/Dynamic/Ida.cpp
+++ b/GridKit/Solver/Dynamic/Ida.cpp
@@ -166,6 +166,10 @@ namespace AnalysisManager
       linearSolver_ = SUNLinSol_KLU(yy_, JacobianMat_, context_);
       checkAllocation((void*) linearSolver_, "SUNLinSol_KLU");
 
+      // set the solver to amd
+      retval = SUNLinSol_KLUSetOrdering(linearSolver_, 0);
+      checkOutput(retval, "SUNLinSol_KLUSetOrdering");
+
       retval = IDASetLinearSolver(solver_, linearSolver_, JacobianMat_);
       checkOutput(retval, "IDASetLinearSolver");
 
diff --git a/GridKit/Utilities/String.hpp b/GridKit/Utilities/String.hpp
index 1a559cf77..d2d402f54 100644
--- a/GridKit/Utilities/String.hpp
+++ b/GridKit/Utilities/String.hpp
@@ -12,14 +12,14 @@ namespace GridKit
     /**
      * @brief Convert a string to all uppercase
      */
-    std::string toUpper(std::string str)
+    inline std::string toUpper(std::string str)
     {
       std::transform(str.begin(), str.end(), str.begin(), [](unsigned char c)
                      { return std::toupper(c); });
       return str;
     }
 
-    std::string strip(std::string str)
+    inline std::string strip(std::string str)
     {
       auto notspace = [](char c)
       { return !std::isspace(c); };