CompFUSE · PDoakORNL · Jun 3, 2026 · May 22, 2026 · May 22, 2026 · May 22, 2026
diff --git a/include/dca/linalg/matrix.hpp b/include/dca/linalg/matrix.hpp
@@ -401,10 +401,15 @@ void Matrix<ScalarType, device_name, ALLOC>::resize(std::pair<int, int> new_size
     // This is slow but still thinking about what to do here, the legacy behavior is fine if
     // we're just overwriting the new sections but not if its assumed they are 0, not all client
     // code is clear about this.
-    if constexpr (device_name == linalg::DeviceType::CPU)
-      for (int i = size_.first; i < new_size.first; ++i)
-        for (int j = size_.second; j < new_size.second; ++j)
-          *(data_ + i * new_size.second + j) = 0;
+    if constexpr (device_name == linalg::DeviceType::CPU) {
+      const int ld = leadingDimension();
+      for (int j = 0; j < new_size.second; ++j)
+        for (int i = size_.first; i < new_size.first; ++i)
+          data_[i + j * ld] = 0;
+      for (int j = size_.second; j < new_size.second; ++j)
+        for (int i = 0; i < size_.first; ++i)
+          data_[i + j * ld] = 0;
+    }
     size_ = new_size;
   }
 #endif

diff --git a/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_walker.hpp b/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_walker.hpp
@@ -518,7 +518,7 @@ CtauxWalker<device_t, Parameters, Data>::CtauxWalker(Parameters& parameters_ref,
   }
 }
 
-  
+
 template <dca::linalg::DeviceType device_t, class Parameters, class Data>
 void CtauxWalker<device_t, Parameters, Data>::printSummary() const {
   // std::defaultfloat is only supported by GCC 5 or later.
@@ -940,11 +940,11 @@ int CtauxWalker<device_t, Parameters, Data>::generateDelayedSpinsAbortAtBennett(
   assert(single_spin_updates_proposed ==
          currently_proposed_creations_ + currently_proposed_annihilations_ + num_statics);
 
-// #ifndef NDEBUG
-//   if (single_spin_updates_proposed >= max_num_delayed_spins) {
-//     std::cout << "single_spin_updates_proposed = " << single_spin_updates_proposed << "  max_num_delayed_spins = " << max_num_delayed_spins << '\n';
-//   }
-// #endif
+#ifndef NDEBUG
+  if (single_spin_updates_proposed > max_num_delayed_spins) {
+    std::cout << "single_spin_updates_proposed = " << single_spin_updates_proposed << "  max_num_delayed_spins = " << max_num_delayed_spins << '\n';
+  }
+#endif
 
   return single_spin_updates_proposed;
 }
@@ -1403,7 +1403,7 @@ void CtauxWalker<device_t, Parameters, Data>::add_delayed_spin(int& delayed_inde
 #ifndef NDEBUG
       isDebugTestMinMax(false, true, ratio_HS_field_DN, Gamma_index_HS_field_DN, Gamma_up_CPU, Gamma_up_diag_max, Gamma_up_diag_min, trace_depth);
 #endif
- 
+
       Gamma_up_size += 1;
     }
     else {
@@ -1459,7 +1459,7 @@ void CtauxWalker<device_t, Parameters, Data>::add_delayed_spin(int& delayed_inde
 	std::cerr << "Gamma_index_HS_field_UP = " << Gamma_index_HS_field_UP << "  Gamma_dn_size = " << Gamma_dn_size << '\n';
       }
 #endif
- 
+
       ratio_HS_field_UP = ctaux_tools.solve_Gamma_blocked(Gamma_index_HS_field_UP, Gamma_dn_CPU,
                                                           exp_delta_V_HS_field_UP,
                                                           Gamma_dn_diag_max, Gamma_dn_diag_min);

diff --git a/include/dca/phys/dca_step/cluster_solver/stdthread_qmci/stdthread_qmci_cluster_solver.hpp b/include/dca/phys/dca_step/cluster_solver/stdthread_qmci/stdthread_qmci_cluster_solver.hpp
@@ -641,7 +641,7 @@ void StdThreadQmciClusterSolver<QmciSolver>::startWalkerAndAccumulator(int id,
   catch (...) {
     throw std::runtime_error("something mysterious  went wrong in walker thread!");
   }
-			     
+
   ++walk_finished_;
   if (BaseClass::writer_ && BaseClass::writer_->isADIOS2())
     BaseClass::writer_->flush();

diff --git a/test/unit/phys/dca_step/cluster_solver/ctaux/structs/CMakeLists.txt b/test/unit/phys/dca_step/cluster_solver/ctaux/structs/CMakeLists.txt
@@ -6,3 +6,10 @@ dca_add_gtest(read_write_config_test
         INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
         LIBS     FFTW::Double ${DCA_LIBS}
         )
+
+dca_add_gtest(ct_aux_hs_configuration_test
+    FAST
+    GTEST_MAIN
+    INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
+    LIBS FFTW::Double ${DCA_LIBS}
+    )
diff --git a/test/unit/phys/dca_step/cluster_solver/ctaux/structs/ct_aux_hs_configuration_test.cpp b/test/unit/phys/dca_step/cluster_solver/ctaux/structs/ct_aux_hs_configuration_test.cpp
@@ -0,0 +1,142 @@
+// Copyright (C) 2025 UT-Battelle, LLC
+// All rights reserved.
+//
+// See LICENSE for terms of usage.
+// See CITATION.md for citation guidelines, if DCA++ is used for scientific publications.
+//
+// Unit tests for CT_AUX_HS_configuration, focusing on get_first_non_interacting_spin_index.
+
+#include "dca/phys/dca_step/cluster_solver/ctaux/structs/ct_aux_hs_configuration.hpp"
+
+#include <vector>
+
+#include "dca/testing/gtest_h_w_warning_blocking.h"
+
+#include "dca/math/random/std_random_wrapper.hpp"
+#include "dca/phys/dca_step/cluster_solver/ctaux/structs/cv.hpp"
+
+using Scalar = double;
+
+#include "test/mock_mcconfig.hpp"
+namespace dca {
+namespace config {
+using McOptions = MockMcOptions<Scalar>;
+}  // namespace config
+}  // namespace dca
+
+#include "test/unit/phys/dca_step/cluster_solver/test_setup.hpp"
+
+constexpr char input_name[] =
+    DCA_SOURCE_DIR "/test/unit/phys/dca_step/cluster_solver/ctaux/structs/input.json";
+
+using CtauxHsConfigurationTest =
+    dca::testing::G0Setup<Scalar, dca::testing::LatticeBilayer, dca::ClusterSolverId::CT_AUX,
+                          input_name>;
+using Parameters = CtauxHsConfigurationTest::Parameters;
+
+TEST_F(CtauxHsConfigurationTest, EmptyConfig) {
+  dca::phys::solver::ctaux::CV<Parameters>::get_H_interaction() = data_->H_interactions;
+
+  std::vector<double> random(200);
+  for (auto& x : random)
+    x = static_cast<double>(std::rand()) / RAND_MAX;
+  Parameters::random_number_generator rng(random);
+  dca::phys::solver::ctaux::CT_AUX_HS_configuration<Parameters> config(parameters_, rng);
+
+  // Empty configuration: method returns 0 (same as size() since both are 0).
+  EXPECT_EQ(0, config.get_first_non_interacting_spin_index(dca::phys::e_UP));
+  EXPECT_EQ(0, config.get_first_non_interacting_spin_index(dca::phys::e_DN));
+}
+
+TEST_F(CtauxHsConfigurationTest, AllAnnihilatableReturnsSize) {
+  dca::phys::solver::ctaux::CV<Parameters>::get_H_interaction() = data_->H_interactions;
+
+  std::vector<double> random(200);
+  for (auto& x : random)
+    x = static_cast<double>(std::rand()) / RAND_MAX;
+  Parameters::random_number_generator rng(random);
+  dca::phys::solver::ctaux::CT_AUX_HS_configuration<Parameters> config(parameters_, rng);
+  config.initialize();
+
+  // After initialize(), all vertices are interacting (annihilatable).
+  // Method should return size() as the "not found" sentinel.
+  const int up_size = config.get(dca::phys::e_UP).size();
+  const int dn_size = config.get(dca::phys::e_DN).size();
+
+  EXPECT_EQ(up_size, config.get_first_non_interacting_spin_index(dca::phys::e_UP));
+  EXPECT_EQ(dn_size, config.get_first_non_interacting_spin_index(dca::phys::e_DN));
+}
+
+TEST_F(CtauxHsConfigurationTest, FirstEntryNonAnnihilatable) {
+  dca::phys::solver::ctaux::CV<Parameters>::get_H_interaction() = data_->H_interactions;
+
+  std::vector<double> random(200);
+  for (auto& x : random)
+    x = static_cast<double>(std::rand()) / RAND_MAX;
+  Parameters::random_number_generator rng(random);
+  dca::phys::solver::ctaux::CT_AUX_HS_configuration<Parameters> config(parameters_, rng);
+  config.initialize();
+
+  // Mark the vertex referenced by configuration_e_UP[0] as non-annihilatable.
+  const auto& up = config.get(dca::phys::e_UP);
+  if (!up.empty()) {
+    int target_idx = up[0].get_configuration_index();
+    config[target_idx].set_annihilatable(false);
+
+    // The first entry is now non-annihilatable, so the method should return 0.
+    EXPECT_EQ(0, config.get_first_non_interacting_spin_index(dca::phys::e_UP));
+    EXPECT_FALSE(config[target_idx].is_annihilatable());
+  }
+
+  // Mark the vertex referenced by configuration_e_DN[0] as non-annihilatable.
+  const auto& dn = config.get(dca::phys::e_DN);
+  if (!dn.empty()) {
+    int target_idx = dn[0].get_configuration_index();
+    config[target_idx].set_annihilatable(false);
+
+    EXPECT_EQ(0, config.get_first_non_interacting_spin_index(dca::phys::e_DN));
+    EXPECT_FALSE(config[target_idx].is_annihilatable());
+  }
+}
+
+TEST_F(CtauxHsConfigurationTest, InsertNoninteractingVertex) {
+  dca::phys::solver::ctaux::CV<Parameters>::get_H_interaction() = data_->H_interactions;
+
+  std::vector<double> random(200);
+  for (auto& x : random)
+    x = static_cast<double>(std::rand()) / RAND_MAX;
+  Parameters::random_number_generator rng(random);
+  dca::phys::solver::ctaux::CT_AUX_HS_configuration<Parameters> config(parameters_, rng);
+  config.initialize();
+
+  int old_size = config.size();
+
+  // Insert a non-interacting vertex (mark_annihilatable = false).
+  config.insert_random_noninteracting_vertex(false);
+
+  EXPECT_EQ(old_size + 1, config.size());
+  EXPECT_FALSE(config[old_size].is_annihilatable());
+
+  // Verify e_UP: the first non-annihilatable index should point to a non-annihilatable vertex,
+  // and all entries before it must be annihilatable.
+  const auto& up = config.get(dca::phys::e_UP);
+  int first_up = config.get_first_non_interacting_spin_index(dca::phys::e_UP);
+  if (first_up < static_cast<int>(up.size())) {
+    int config_idx = up[first_up].get_configuration_index();
+    EXPECT_FALSE(config[config_idx].is_annihilatable());
+    for (int i = 0; i < first_up; ++i) {
+      EXPECT_TRUE(config[up[i].get_configuration_index()].is_annihilatable());
+    }
+  }
+
+  // Same check for e_DN.
+  const auto& dn = config.get(dca::phys::e_DN);
+  int first_dn = config.get_first_non_interacting_spin_index(dca::phys::e_DN);
+  if (first_dn < static_cast<int>(dn.size())) {
+    int config_idx = dn[first_dn].get_configuration_index();
+    EXPECT_FALSE(config[config_idx].is_annihilatable());
+    for (int i = 0; i < first_dn; ++i) {
+      EXPECT_TRUE(config[dn[i].get_configuration_index()].is_annihilatable());
+    }
+  }
+}
diff --git a/tutorials/tc/gen_temps.awk b/tutorials/tc/gen_temps.awk
@@ -0,0 +1,94 @@
+#!/usr/bin/awk -f
+# Generates the tc/preconfigured input files from the templates input_sp.json.in and input_tp.json.in.
+# Run this script from tutorials/tc/.
+
+BEGIN {
+    # Temperature series from high to low.
+    n_temps = 9
+    temps[1]  = 1
+    temps[2]  = 0.75
+    temps[3]  = 0.5
+    temps[4]  = 0.25
+    temps[5]  = 0.125
+    temps[6]  = 0.1
+    temps[7]  = 0.09
+    temps[8]  = 0.08
+    temps[9]  = 0.07
+
+    # Pre-computed beta = 1/T values matching the existing preconfigured files.
+    betas[1]  = 1
+    betas[2]  = 1.33333
+    betas[3]  = 2
+    betas[4]  = 4
+    betas[5]  = 8
+    betas[6]  = 10
+    betas[7]  = 11.1111
+    betas[8]  = 12.5
+    betas[9]  = 14.2857
+
+    # Temperatures that also need a two-particle (tp) input file.
+    has_tp[6] = 1
+    has_tp[7] = 1
+    has_tp[8] = 1
+    has_tp[9] = 1
+
+    # Fixed physical parameters.
+    dens = 0.95
+    hubbardU = 6
+    vec1 = "[2, 0]"
+    vec2 = "[0, 2]"
+
+    for (i = 1; i <= n_temps; i++) {
+        t = temps[i]
+
+        # Previous temperature for self-energy restart ("zero" for the highest T).
+        prev = (i == 1) ? "zero" : temps[i - 1]
+
+        # DCA iterations: 8 for the first temperature, 6 otherwise.
+        iters_sp = (i == 1) ? 8 : 6
+
+        dir = "preconfigured/T=" t
+        cmd = "mkdir -p " dir
+        system(cmd)
+
+        # Generate input_sp.json.
+        sp_in = "input_sp.json.in"
+        sp_out = dir "/input_sp.json"
+        while ((getline line < sp_in) > 0) {
+            gsub(/BETA/, betas[i], line)
+            gsub(/DENS/, dens, line)
+            gsub(/HUBBARDU/, hubbardU, line)
+            gsub(/VEC1/, vec1, line)
+            gsub(/VEC2/, vec2, line)
+            gsub(/CURRENT_TEMP/, t, line)
+            gsub(/PREVIOUS_TEMP/, prev, line)
+            # Fix the initial self-energy for the highest temperature.
+            gsub(/"\.\/T=zero\/dca_sp\.hdf5"/, "\"zero\"", line)
+            gsub(/ITERS/, iters_sp, line)
+            print line > sp_out
+        }
+        close(sp_in)
+        close(sp_out)
+
+        # Generate input_tp.json for temperatures that need it.
+        if (has_tp[i]) {
+            tp_in = "input_tp.json.in"
+            tp_out = dir "/input_tp.json"
+            while ((getline line < tp_in) > 0) {
+                gsub(/BETA/, betas[i], line)
+                gsub(/DENS/, dens, line)
+                gsub(/HUBBARDU/, hubbardU, line)
+                gsub(/VEC1/, vec1, line)
+                gsub(/VEC2/, vec2, line)
+                gsub(/CURRENT_TEMP/, t, line)
+                # PREVIOUS_TEMP is not used in the tp template, but substitute
+                # with the current temperature for safety.
+                gsub(/PREVIOUS_TEMP/, t, line)
+                gsub(/ITERS/, 1, line)
+                print line > tp_out
+            }
+            close(tp_in)
+            close(tp_out)
+        }
+    }
+}
diff --git a/tutorials/tc/input_sp.json.in b/tutorials/tc/input_sp.json.in
@@ -62,16 +62,17 @@
         "measurements": 100000,
 
         "threaded-solver": {
-            "walkers": 3,
-            "accumulators": 5
+            "walkers": 4,
+            "accumulators": 4,
+            "shared-walk-and-accumulation-thread": true
         }
     },
 
     "CT-AUX": {
         "expansion-parameter-K": 1.,
         "initial-configuration-size": 16,
         "initial-matrix-size": 16,
-        "max-submatrix-size": 16,
+        "max-submatrix-size": 256,
         "neglect-Bennett-updates": false,
         "additional-time-measurements": false
     }

diff --git a/tutorials/tc/input_tp.json.in b/tutorials/tc/input_tp.json.in
@@ -64,16 +64,17 @@
         "measurements": 100000,
 
         "threaded-solver": {
-            "walkers": 3,
-            "accumulators": 5
+            "walkers": 4,
+            "accumulators": 4,
+            "shared-walk-and-accumulation-thread": true
         }
     },
 
     "CT-AUX": {
         "expansion-parameter-K": 1.,
         "initial-configuration-size": 16,
         "initial-matrix-size": 16,
-        "max-submatrix-size": 16,
+        "max-submatrix-size": 256,
         "neglect-Bennett-updates": false,
         "additional-time-measurements": false
     },

diff --git a/tutorials/tc/preconfigured/README.md b/tutorials/tc/preconfigured/README.md
@@ -1,5 +1,5 @@
 This directory contains ready-to-use input files and jobs scripts for the [T<sub>c</sub> tutorial](https://github.com/CompFUSE/DCA/wiki/Tutorial:-Tc).
-The files are preconfigured to run the example with a single process (no MPI) using the threaded Monte Carlo solver with 3 walkers and 5 accumulators.
+The files are preconfigured to run the example with a single process (no MPI) using the threaded Monte Carlo solver with 1 walker and 1 accumulator per MPI rank.
 The total number of measurements is 100 000.
 
 The temperature steps are given by

diff --git a/tutorials/tc/preconfigured/T=0.07/input_sp.json b/tutorials/tc/preconfigured/T=0.07/input_sp.json
@@ -62,16 +62,17 @@
         "measurements": 100000,
 
         "threaded-solver": {
-            "walkers": 3,
-            "accumulators": 5
+            "walkers": 4,
+            "accumulators": 4,
+            "shared-walk-and-accumulation-thread": true
         }
     },
 
     "CT-AUX": {
         "expansion-parameter-K": 1.,
         "initial-configuration-size": 16,
         "initial-matrix-size": 16,
-        "max-submatrix-size": 16,
+        "max-submatrix-size": 256,
         "neglect-Bennett-updates": false,
         "additional-time-measurements": false
     }