Refactor ctint walker

cmake/dca_config.cmake

@@ -145,6 +145,8 @@ elseif (DCA_LATTICE STREQUAL "Kagome")
     "dca/phys/models/analytic_hamiltonians/Kagome_hubbard.hpp")
 elseif (DCA_LATTICE STREQUAL "hund")
   set(DCA_LATTICE_TYPE dca::phys::models::HundLattice<PointGroup>)
+  set(DCA_LATTICE_INCLUDE
+    "dca/phys/models/analytic_hamiltonians/hund_lattice.hpp")
 elseif (DCA_LATTICE STREQUAL "threeband")
   set(DCA_LATTICE_TYPE dca::phys::models::ThreebandHubbard<PointGroup>)
   set(DCA_LATTICE_INCLUDE

include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_cluster_solver.hpp

@@ -117,6 +117,8 @@ class CtauxClusterSolver {
     return g0_;
   };
 
+  Walker::Resource& getResource() { return dummy_walker_resource_; };
+
 protected:
   void warmUp(Walker& walker);
 
@@ -161,6 +163,7 @@ class CtauxClusterSolver {
 
   G0Interpolation<device, Scalar> g0_;
 
+  Walker::Resource walker_dummy_resource_;
 private:
   Rng rng_;
 

include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_walker.hpp

@@ -66,6 +66,9 @@ class CtauxWalker : public WalkerBIT<Parameters, Data>, public CtauxWalkerData<d
   constexpr static dca::linalg::DeviceType device = device_t;
   static constexpr bool is_complex = dca::util::IsComplex<Scalar>::value;
 
+  struct DummyResource {};
+  using Resource = DummyResource;
+
 public:
   /** constructor
    *  param[in] id    thread id

include/dca/phys/dca_step/cluster_solver/ctint/ctint_cluster_solver.hpp

@@ -35,6 +35,7 @@
 #include "dca/phys/dca_step/cluster_solver/ctint/details/solver_methods.hpp"
 #include "dca/phys/dca_step/cluster_solver/ctint/domains/common_domains.hpp"
 #include "dca/phys/dca_step/cluster_solver/ctint/walker/ctint_walker_choice.hpp"
+#include "dca/phys/dca_step/cluster_solver/ctint/walker/tools/d_matrix_builder.hpp"
 #include "dca/phys/dca_step/cluster_solver/shared_tools/interpolation/g0_interpolation.hpp"
 //#include "dca/phys/dca_step/cluster_solver/shared_tools/accumulation/time_correlator.hpp"
 #include "dca/phys/dca_data/dca_data.hpp"
@@ -74,6 +75,8 @@ class CtintClusterSolver {
   using Data = DcaData<Parameters, DIST>;
   static constexpr linalg::DeviceType device = device_t;
 
+  using DMatrixBuilder = ctint::DMatrixBuilder<device_t, Scalar>;
+
   CtintClusterSolver(Parameters& parameters_ref, Data& Data_ref,
                      std::shared_ptr<io::Writer<Concurrency>> writer);
 
@@ -108,7 +111,8 @@ class CtintClusterSolver {
   // Returns the function G(k,w) without averaging across MPI ranks.
   auto local_G_k_w() const;
 
-protected:  // thread jacket interface.
+  DMatrixBuilder& getResource() { return *d_matrix_builder_; };
+protected:  // thread jacket interface.  
   using ParametersType = Parameters;
   using DataType = Data;
   using Rng = typename Parameters::random_number_generator;
@@ -166,6 +170,7 @@ class CtintClusterSolver {
   std::unique_ptr<Walker> walker_;
   // Walker input.
   Rng rng_;
+  std::unique_ptr<DMatrixBuilder> d_matrix_builder_;
 };
 
 template <dca::linalg::DeviceType DEV, class PARAM, bool use_submatrix, DistType DIST>
@@ -174,11 +179,10 @@ CtintClusterSolver<DEV, PARAM, use_submatrix, DIST>::CtintClusterSolver(
     : parameters_(parameters_ref),
       concurrency_(parameters_.get_concurrency()),
       data_(data_ref),
-
       accumulator_(parameters_, data_),
       writer_(writer),
       rng_(concurrency_.id(), concurrency_.number_of_processors(), parameters_.get_seed()) {
-  Walker::setDMatrixBuilder(g0_);
+  d_matrix_builder_ = std::make_unique<DMatrixBuilder>(g0_, PARAM::bands, RDmn());
   Walker::setInteractionVertices(data_, parameters_);
 
   if (concurrency_.id() == concurrency_.first())
@@ -191,7 +195,7 @@ void CtintClusterSolver<device_t, Parameters, use_submatrix, DIST>::initialize(i
 
   g0_.initializeShrinked(data_.G0_r_t_cluster_excluded);
 
-  Walker::setDMatrixAlpha(parameters_.getAlphas(), parameters_.adjustAlphaDd());
+  d_matrix_builder_->setAlphas(parameters_.getAlphas(), parameters_.adjustAlphaDd());
 
   // It is a waiting to happen bug for this to be here and in CtintAccumulator
   accumulator_.initialize(dca_iteration_);
@@ -203,7 +207,7 @@ void CtintClusterSolver<device_t, Parameters, use_submatrix, DIST>::initialize(i
 
 template <dca::linalg::DeviceType device_t, class Parameters, bool use_submatrix, DistType DIST>
 void CtintClusterSolver<device_t, Parameters, use_submatrix, DIST>::integrate() {
-  walker_ = std::make_unique<Walker>(parameters_, data_, rng_, 0);
+  walker_ = std::make_unique<Walker>(parameters_, data_, rng_, d_matrix_builder_, 0);
   walker_->initialize(dca_iteration_);
 
   dca::profiling::WallTime start_time;

include/dca/phys/dca_step/cluster_solver/ctint/walker/ctint_walker_base.hpp

@@ -80,15 +80,22 @@ class CtintWalkerBase {
   void computeM(MatrixPair& m_accum) const;
 
   // Reset the counters and recompute the configuration sign and weight.
-  void markThermalized();
+  virtual void markThermalized() = 0;
 
   /** Recompute the matrix M from the configuration in O(expansion_order^3) time.
    *  Postcondition: sign_ and mc_log_weight_ are recomputed.
    *  mc_log_weight is the negative sum of the log det of both sectors of M
    *  + log of the abs of each vertices interaction strength.
    */
-  virtual void setMFromConfig();
+  virtual void setMFromConfig() = 0;
 
+  virtual void doStep() = 0;
+
+  virtual void doSweep() = 0;
+
+  template<linalg::DeviceType DEVICE>
+  void setMFromConfigImpl(DMatrixBuilder<DEVICE, Scalar>& d_matrix_builder);
+
   bool is_thermalized() const {
     return thermalized_;
   }
@@ -152,14 +159,6 @@ class CtintWalkerBase {
     buff >> configuration_;
   }
 
-  // Initialize the builder object shared by all walkers.
-  template <linalg::DeviceType device_type>
-  static void setDMatrixBuilder(const G0Interpolation<device_type, Scalar>& g0);
-
-  static void setDMatrixAlpha(const std::array<double, 3>& alphas, bool adjust_dd);
-
-  static void setInteractionVertices(const Data& data, const Parameters& parameters);
-
   double stealFLOPs() {
     auto flop = flop_;
     flop_ = 0.;
@@ -173,12 +172,11 @@ class CtintWalkerBase {
 
   static void sumConcurrency(const Concurrency&) {}
 
-  const auto& get_dmatrix_builder() const {
-    return *d_builder_ptr_;
-  }
-
   void writeAlphas() const;
 
+  static void setInteractionVertices(const Data& data,
+			      const Parameters& parameters);
+
 protected:
   // typedefs
   using RDmn = typename Parameters::RClusterDmn;
@@ -187,7 +185,6 @@ class CtintWalkerBase {
   void updateSweepAverages();
 
 protected:  // Members.
-  static inline std::unique_ptr<DMatrixBuilder<linalg::CPU, Scalar>> d_builder_ptr_;
   static inline InteractionVertices vertices_;
 
   const Parameters& parameters_;
@@ -271,41 +268,6 @@ void CtintWalkerBase<Parameters,DIST>::initialize(int iteration) {
   setMFromConfig();
 }
 
-template <class Parameters, DistType DIST>
-void CtintWalkerBase<Parameters,DIST>::setMFromConfig() {
-  mc_log_weight_ = 0.;
-  phase_.reset();
-
-  for (int s = 0; s < 2; ++s) {
-    // compute Mij = g0(t_i,t_j) - I* alpha(s_i)
-
-    const auto& sector = configuration_.getSector(s);
-    auto& M = M_[s];
-    const int n = sector.size();
-    M.resize(n);
-    if (!n)
-      continue;
-    for (int j = 0; j < n; ++j)
-      for (int i = 0; i < n; ++i)
-        M(i, j) = d_builder_ptr_->computeD(i, j, sector);
-
-    if (M.nrRows()) {
-      const auto [log_det, phase] = linalg::matrixop::inverseAndLogDeterminant(M);
-
-      mc_log_weight_ -= log_det;  // Weight proportional to det(M^{-1})
-      phase_.divide(phase);
-    }
-  }
-
-  // So what is going on here.
-  for (int i = 0; i < configuration_.size(); ++i) {
-    // This is actual interaction strength of the vertex i.e H_int(nu1, nu2, delta_r)
-    const Real term = -configuration_.getStrength(i);
-    mc_log_weight_ += std::log(std::abs(term));
-    phase_.multiply(term);
-  }
-}
-
 template <class Parameters, DistType DIST>
 void CtintWalkerBase<Parameters,DIST>::updateSweepAverages() {
   order_avg_.addSample(order());
@@ -315,39 +277,22 @@ void CtintWalkerBase<Parameters,DIST>::updateSweepAverages() {
     partial_order_avg_.addSample(order());
 }
 
-template <class Parameters, DistType DIST>
-void CtintWalkerBase<Parameters,DIST>::writeAlphas() const {
-  std::cout << "For initial configuration integration:\n";
-  for (int isec = 0; isec < 2; ++isec) {
-    std::cout << "Sector: " << isec << '\n';
-    for (int ic = 0; ic < order(); ++ic) {
-      auto aux_spin = configuration_.getSector(isec).getAuxFieldType(ic);
-      auto left_b = configuration_.getSector(isec).getLeftB(ic);
-      auto alpha_left = get_dmatrix_builder().computeAlpha(aux_spin, left_b);
-      std::cout << "vertex: " << std::setw(6) << ic;
-      std::cout << " | aux spin: " << aux_spin << " | left B: " << left_b
-                << " | alpha left = " << alpha_left << '\n';
-    }
-  }
-}
-
-template <class Parameters, DistType DIST>
-void CtintWalkerBase<Parameters,DIST>::markThermalized() {
-  thermalized_ = true;
-
-  nb_steps_per_sweep_ = std::max(1., std::ceil(sweeps_per_meas_ * partial_order_avg_.mean()));
-  thermalization_steps_ = n_steps_;
-
-  order_avg_.reset();
-  sign_avg_.reset();
-  n_accepted_ = 0;
+// template <class Parameters, DistType DIST>
+// void CtintWalkerBase<Parameters,DIST>::writeAlphas() const {
+//   std::cout << "For initial configuration integration:\n";
+//   for (int isec = 0; isec < 2; ++isec) {
+//     std::cout << "Sector: " << isec << '\n';
+//     for (int ic = 0; ic < order(); ++ic) {
+//       auto aux_spin = configuration_.getSector(isec).getAuxFieldType(ic);
+//       auto left_b = configuration_.getSector(isec).getLeftB(ic);
+//       auto alpha_left = get_dmatrix_builder().computeAlpha(aux_spin, left_b);
+//       std::cout << "vertex: " << std::setw(6) << ic;
+//       std::cout << " | aux spin: " << aux_spin << " | left B: " << left_b
+//                 << " | alpha left = " << alpha_left << '\n';
+//     }
+//   }
+// }
 
-  // Recompute the Monte Carlo weight.
-  setMFromConfig();
-#ifndef NDEBUG
-  //writeAlphas();
-#endif
-}
 
 template <class Parameters, DistType DIST>
 void CtintWalkerBase<Parameters,DIST>::updateShell(int meas_id, int meas_to_do) const {
@@ -383,25 +328,6 @@ void CtintWalkerBase<Parameters,DIST>::printSummary() const {
   std::cout << std::endl;
 }
 
-template <class Parameters, DistType DIST>
-template <linalg::DeviceType device_type>
-void CtintWalkerBase<Parameters,DIST>::setDMatrixBuilder(
-							 const dca::phys::solver::G0Interpolation<device_type, Scalar>& g0) {
-  using RDmn = typename Parameters::RClusterDmn;
-
-  if (d_builder_ptr_)
-    std::cerr << "Warning: DMatrixBuilder already set." << std::endl;
-
-  d_builder_ptr_ = std::make_unique<DMatrixBuilder<device_type, Scalar>>(g0, n_bands_, RDmn());
-}
-
-template <class Parameters, DistType DIST>
-void CtintWalkerBase<Parameters,DIST>::setDMatrixAlpha(const std::array<double, 3>& alphas,
-                                                              bool adjust_dd) {
-  assert(d_builder_ptr_);
-  d_builder_ptr_->setAlphas(alphas, adjust_dd);
-}
-
 template <class Parameters, DistType DIST>
 void CtintWalkerBase<Parameters,DIST>::setInteractionVertices(const Data& data,
                                                                      const Parameters& parameters) {
@@ -419,6 +345,78 @@ void CtintWalkerBase<Parameters,DIST>::computeM(MatrixPair& m_accum) const {
   m_accum = M_;
 }
 
+// template<class WALKER, linalg::DeviceType DEVICE>  
+// void setMFromConfigHelper(WALKER& walker, DMatrixBuilder<DEVICE, Scalar>& d_matrix_builder) {
+//   walker.mc_log_weight_ = 0.;
+//   walker.phase_.reset();
+
+//   for (int s = 0; s < 2; ++s) {
+//     // compute Mij = g0(t_i,t_j) - I* alpha(s_i)
+
+//     const auto& sector = walker.configuration_.getSector(s);
+//     auto& M = walker.M_[s];
+//     const int n = sector.size();
+//     M.resize(n);
+//     if (!n)
+//       continue;
+//     for (int j = 0; j < n; ++j)
+//       for (int i = 0; i < n; ++i)
+//         M(i, j) = d_matrix_builder.computeD(i, j, sector);
+
+//     if (M.nrRows()) {
+//       const auto [log_det, phase] = linalg::matrixop::inverseAndLogDeterminant(M);
+
+//       walker.mc_log_weight_ -= log_det;  // Weight proportional to det(M^{-1})
+//       walker.phase_.divide(phase);
+//     }
+//   }
+
+//   // So what is going on here.
+//   for (int i = 0; i < walker.configuration_.size(); ++i) {
+//     // This is actual interaction strength of the vertex i.e H_int(nu1, nu2, delta_r)
+//     const typename decltype(walker)::Real term = -walker.configuration_.getStrength(i);
+//     walker.mc_log_weight_ += std::log(std::abs(term));
+//     walker.phase_.multiply(term);
+//   }
+// }
+
+
+template <class Parameters, DistType DIST>
+template <linalg::DeviceType DEVICE>
+void CtintWalkerBase<Parameters, DIST>::setMFromConfigImpl(DMatrixBuilder<DEVICE, Scalar>& d_matrix_builder) {
+  mc_log_weight_ = 0.;
+  phase_.reset();
+
+  for (int s = 0; s < 2; ++s) {
+    // compute Mij = g0(t_i,t_j) - I* alpha(s_i)
+
+    const auto& sector = configuration_.getSector(s);
+    auto& M = M_[s];
+    const int n = sector.size();
+    M.resize(n);
+    if (!n)
+      continue;
+    for (int j = 0; j < n; ++j)
+      for (int i = 0; i < n; ++i)
+        M(i, j) = d_matrix_builder.computeD(i, j, sector);
+
+    if (M.nrRows()) {
+      const auto [log_det, phase] = linalg::matrixop::inverseAndLogDeterminant(M);
+
+      mc_log_weight_ -= log_det;  // Weight proportional to det(M^{-1})
+      phase_.divide(phase);
+    }
+  }
+
+  // So what is going on here.
+  for (int i = 0; i < configuration_.size(); ++i) {
+    // This is actual interaction strength of the vertex i.e H_int(nu1, nu2, delta_r)
+    const Real term = -configuration_.getStrength(i);
+    mc_log_weight_ += std::log(std::abs(term));
+    phase_.multiply(term);
+  }
+}
+
 }  // namespace ctint
 }  // namespace solver
 }  // namespace phys