It now compiles in a non-MPI environment (not well tested)

Author: shinaoka

CMakeLists.txt

@@ -8,7 +8,9 @@ list(APPEND CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/cmake ${CMAKE_SOURCE_DIR})
 # export ALPSCore_DIR=/location/to/ALPSCORE/
 find_package(ALPSCore REQUIRED)
 find_package(Boost 1.55.0 REQUIRED timer)
-find_package(MPI REQUIRED)
+if (ALPS_ENABLE_MPI)
+    find_package(MPI REQUIRED)
+endif()
 find_package(Eigen3 REQUIRED)
 
 # Option (use quad precision for part of calculations)

src/impurity.hpp

@@ -39,7 +39,9 @@
 #include <alps/accumulators.hpp>
 #include <alps/mc/api.hpp>
 #include <alps/mc/mcbase.hpp>
-#include <alps/mc/mpiadapter.hpp>
+#ifdef ALPS_HAVE_MPI
+  #include <alps/mc/mpiadapter.hpp>
+#endif
 #include <alps/mc/stop_callback.hpp>
 #include <alps/params/convenience_params.hpp>
 
@@ -108,7 +110,9 @@ class HybridizationSimulation : public alps::mcbase
   boost::scoped_ptr<IMP_MODEL> p_model;
 
   //ALPS MPI communicator
+#ifdef ALPS_HAVE_MPI
   alps::mpi::communicator comm;
+#endif
 
   //Constant simulation parameters
   const long thermalization_sweeps;           // sweeps to be done for equilibration

src/impurity.ipp

@@ -48,7 +48,9 @@ HybridizationSimulation<IMP_MODEL>::HybridizationSimulation(parameters_type cons
     N(static_cast<int>(parameters["N_TAU"])),                  //time slices
     Np1(N+1),
     p_model(new IMP_MODEL(p,rank==0)),//impurity model
+#ifdef ALPS_HAVE_MPI
     comm(),
+#endif
     thermalization_sweeps(parameters["THERMALIZATION"]),          //sweeps needed for thermalization
     total_sweeps(parameters["SWEEPS"]),                           //sweeps needed for total run
     N_meas(parameters["N_MEAS"]),
@@ -513,18 +515,24 @@ void HybridizationSimulation<IMP_MODEL>::update_MC_parameters() {
     1.0/std::max(static_cast<double>(sweeps)/static_cast<double>(interval_update_cutoff),1.0)
   );
 
+#ifdef ALPS_HAVE_MPI
   boost::tuple<bool,double> r_pair = weight_vs_distance.update_cutoff(acc_rate_cutoff, max_distance_pair, mag, comm);
+#else
+  boost::tuple<bool,double> r_pair = weight_vs_distance.update_cutoff(acc_rate_cutoff, max_distance_pair, mag);
+#endif
   max_distance_pair = boost::get<1>(r_pair);
   max_distance_pair = std::min(0.5*BETA, max_distance_pair);
 
   //std::cout << "Done max_distance_pair rank " << comm.rank() << " sweeps " << sweeps << std::endl;
 
+#ifdef ALPS_HAVE_MPI
   boost::tuple<bool,double> r_shift = weight_vs_distance_shift.update_cutoff(acc_rate_cutoff, max_distance_shift, mag, comm);
+#else
+  boost::tuple<bool,double> r_shift = weight_vs_distance_shift.update_cutoff(acc_rate_cutoff, max_distance_shift, mag);
+#endif
   max_distance_shift = boost::get<1>(r_shift);
   max_distance_shift = std::min(0.5*BETA, max_distance_shift);
 
-  //std::cout << "Done max_distance_shift rank " << comm.rank() << " sweeps " << sweeps << std::endl;
-
   const double max_distance = std::max(max_distance_pair,max_distance_shift);
   const std::size_t n_window_new = static_cast<std::size_t>(std::max(1,static_cast<int>(BETA/(2.0*max_distance))));
 
@@ -544,16 +552,15 @@ void HybridizationSimulation<IMP_MODEL>::prepare_for_measurement() {
   max_dist_optimizer.reset();
   weight_vs_distance.reset();
   weight_vs_distance_shift.reset();
-  //std::cout << "Call prepare_for_measurement rank " << comm.rank() << std::endl;
-  if (comm.rank()==0) {
+  if (global_mpi_rank==0) {
     std::cout << "We're done with thermalization." << std::endl << "The number of segments for sliding window update is " << sliding_window.get_n_window() << "." << std::endl << std::endl;
   }
 
   //N_meas
   const int N_meas_min = std::max(10, 4*sliding_window.get_n_window());//a sweep of the window takes 4*get_n_window()
   if (N_meas<N_meas_min) {
     N_meas = N_meas_min;
-    if (comm.rank()==0) {
+    if (global_mpi_rank==0) {
       std::cout << "Warning N_MEAS is too small: using N_MEAS = "  << N_meas << " instead." << std::endl;
     }
   }
@@ -562,7 +569,7 @@ void HybridizationSimulation<IMP_MODEL>::prepare_for_measurement() {
   const int N_meas_g_min = std::max(10, sliding_window.get_n_window());//a sweep of the window takes 4*get_n_window()
   if (N_meas_g<N_meas_g_min) {
     N_meas_g = N_meas_g_min;
-    if (comm.rank()==0) {
+    if (global_mpi_rank==0) {
       std::cout << "Warning N_MEAS_GREENS_FUNCTION is too small: using N_MEAS_GREENS_FUNCTION = "  << N_meas_g << " instead." << std::endl;
     }
   }
@@ -571,7 +578,7 @@ void HybridizationSimulation<IMP_MODEL>::prepare_for_measurement() {
   const int N_swap_min =std::max(10, 4*sliding_window.get_n_window());//a sweep of the window takes 4*get_n_window()
   if (N_swap<N_swap_min && swap_vector.size()>0) {
     N_swap = N_swap_min;
-    if (comm.rank()==0) {
+    if (global_mpi_rank==0) {
       std::cout << "Warning N_SWAP is too small: using N_SWAP = "  << N_swap << " instead." << std::endl;
     }
   }

src/impurity_init.ipp

@@ -97,7 +97,7 @@ void HybridizationSimulation<IMP_MODEL>::resize_vectors() {
     }
     swap_acc_rate.resize(swap_vector.size());
 
-    if (comm.rank() == 0) {
+    if (global_mpi_rank == 0) {
       std::cout << "The following swap updates will be performed." << std::endl;
       for (int i=0; i<swap_vector.size(); ++i) {
         std::cout << "Update #" << i << " generated from template #" << swap_vector[i].second << std::endl;
@@ -125,7 +125,7 @@ void HybridizationSimulation<IMP_MODEL>::read_eq_time_two_particle_greens_meas()
   const int GF_RANK=2;
 
   const std::string fname_key = "EQUAL_TIME_TWO_PARTICLE_GREENS_FUNCTION";
-  const bool verbose = (comm.rank()==0);
+  const bool verbose = (global_mpi_rank==0);
 
   if (!par.defined(fname_key)) {
     return;
@@ -170,7 +170,7 @@ void HybridizationSimulation<IMP_MODEL>::read_eq_time_two_particle_greens_meas()
 template<typename IMP_MODEL>
 void HybridizationSimulation<IMP_MODEL>::read_two_time_correlation_functions() {
   const std::string fname_key = "TWO_TIME_CORRELATION_FUNCTIONS";
-  const bool verbose = (comm.rank()==0);
+  const bool verbose = (global_mpi_rank==0);
 
   if (!par.defined("N_TAU_TWO_TIME_CORRELATION_FUNCTIONS")) {
     return;

src/main.hpp

@@ -16,51 +16,71 @@
 #include "util.hpp"
 
 #include <alps/utilities/signal.hpp>
-#include <alps/utilities/mpi.hpp>
 #include <alps/mc/api.hpp>
 #include <alps/mc/mcbase.hpp>
-#include <alps/mc/mpiadapter.hpp>
 #include <alps/mc/stop_callback.hpp>
-
-#include "mympiadapter.hpp"
+#ifdef ALPS_HAVE_MPI
+  #include <alps/utilities/mpi.hpp>
+  #include "mc/mympiadapter.hpp"
+#endif
+  #include "mc/mymcadapter.hpp"
 #include "postprocess.hpp"
 
 #undef BUILD_PYTHON_MODULE
 
+int global_mpi_rank;
+
 template<class SOLVER_TYPE>
 int run_simulation(int argc, const char* argv[], typename alps::parameters_type<SOLVER_TYPE>::type& parameters) {
+#ifdef ALPS_HAVE_MPI
   typedef mymcmpiadapter<SOLVER_TYPE> sim_type;
+#else
+  typedef mymcadapter<SOLVER_TYPE> sim_type;
+#endif
 
   SOLVER_TYPE::define_parameters(parameters);
   if (parameters.help_requested(std::cout)) {
     exit(0);
   }
 
   char** argv_tmp = const_cast<char**>(argv);//ugly solution
+#ifdef ALPS_HAVE_MPI
   alps::mpi::environment env(argc, argv_tmp);
   alps::mpi::communicator c;
   c.barrier();
-  if (c.rank()==0) {
+  global_mpi_rank = c.rank();
+  if (global_mpi_rank == 0) {
     std::cout << "Creating simulation..." << std::endl;
   }
   sim_type sim(parameters, c);
-
-  // Run the simulation
   const boost::function<bool()> cb = alps::stop_callback(c,size_t(parameters["timelimit"]));
+#else
+  global_mpi_rank = 0;
+  sim_type sim(parameters);
+  const boost::function<bool()> cb = alps::stop_callback(size_t(parameters["timelimit"]));
+#endif
+
   sim.run(cb);
 
   // Saving to the output file
+#ifdef ALPS_HAVE_MPI
   if (c.rank()==0){
+#endif
     typename alps::results_type<SOLVER_TYPE>::type results = alps::collect_results(sim);
     std::string output_file = parameters["outputfile"];
     alps::hdf5::archive ar(boost::filesystem::path(output_file), "w");
     ar["/parameters"] << parameters;
     ar["/simulation/results"] << results;
     compute_greens_functions<SOLVER_TYPE>(results, parameters, ar);
+#ifdef ALPS_HAVE_MPI
   } else{
     alps::collect_results(sim);
   }
+#endif
+
+#ifdef ALPS_HAVE_MPI
   c.barrier();
+#endif
 
   return 0;
 }

src/mc/mymcadapter.hpp

@@ -0,0 +1,57 @@
+#pragma once
+
+#include <boost/format.hpp>
+#include <boost/date_time/posix_time/posix_time_types.hpp>
+#include <alps/mc/check_schedule.hpp>
+
+template<typename Base, typename ScheduleChecker = alps::check_schedule> class mymcadapter : public Base {
+public:
+  typedef typename Base::parameters_type parameters_type;
+  typedef boost::posix_time::ptime ptime;
+
+  /// Construct mcmpiadapter with alps::check_schedule with the relevant parameters Tmin and Tmax taken from the provided parameters
+  mymcadapter(parameters_type const & parameters)
+    : Base(parameters, 0),
+      schedule_checker(parameters["Tmin"], parameters["Tmax"]) {}
+
+  bool run(boost::function<bool()> const &stop_callback) {
+    bool done = false, stopped = false;
+    bool was_thermalized_before = false;
+    const ptime start_time = boost::posix_time::second_clock::local_time();
+    ptime time_last_output = start_time;
+    const double min_output_interval = 1; //1 sec
+    do {
+      const bool is_thermalized = this->is_thermalized();
+      if (is_thermalized && !was_thermalized_before) {
+        this->prepare_for_measurement();
+      }
+
+      this->update();
+
+      if (is_thermalized) {
+        this->measure();
+      }
+
+      was_thermalized_before = is_thermalized;
+      const ptime current_time = boost::posix_time::second_clock::local_time();
+      if (!stopped && !is_thermalized) {
+        if ((current_time-time_last_output).total_seconds()>min_output_interval) {
+          std::cout << boost::format("Not thermalized yet: %1% sec passed.")%static_cast<int>((current_time-start_time).total_seconds()) << std::endl;
+          time_last_output = current_time;
+        }
+      } else if (stopped || schedule_checker.pending()) {
+        stopped = stop_callback();
+        double local_fraction = stopped ? 1. : Base::fraction_completed();
+        schedule_checker.update(local_fraction);
+        if ((current_time-time_last_output).total_seconds()>min_output_interval) {
+          std::cout << "Checking if the simulation is finished: " << std::min(static_cast<int>(local_fraction*100),100) << "% of Monte Carlo steps done." << std::endl;
+          time_last_output = current_time;
+        }
+        done = local_fraction >= 1.;
+      }
+    } while (!done);
+    return !stopped;
+  }
+private:
+  ScheduleChecker schedule_checker;
+};

src/mympiadapter.hpp src/mc/mympiadapter.hppsrc/mympiadapter.hpp renamed to src/mc/mympiadapter.hpp

@@ -2,6 +2,7 @@
 
 #include <algorithm>
 
+#ifdef ALPS_HAVE_MPI
 #include <alps/utilities/mpi.hpp>
 
 /// performs mpi_allreduce() for a valarray of type T.
@@ -12,6 +13,7 @@ void all_reduce(const alps::mpi::communicator& comm, const std::valarray<T>& in_
   out_vals.resize(in_vals.size());
   MPI_Allreduce((void*)&in_vals[0], (void*)&out_vals[0], in_vals.size(), alps::mpi::detail::mpi_type<T>(), alps::mpi::is_mpi_op<Op,T>::op(), comm);
 }
+#endif
 
 template<class T>
 void rebin(std::valarray<T>& org_array, int nrebin) {
@@ -106,7 +108,11 @@ class scalar_histogram
   }
 
   //maxdist is not updated if we do not have enough data.
+#ifdef ALPS_HAVE_MPI
   boost::tuple<bool,double> update_cutoff(double cutoff_ratio, double maxdist, double mag, const alps::mpi::communicator& alps_comm) const {
+#else
+  boost::tuple<bool,double> update_cutoff(double cutoff_ratio, double maxdist, double mag) const {
+#endif
     assert(cutoff_ratio>=0.0 && cutoff_ratio<=1.0);
     assert(mag>=1.0);
     const int min_count = 10;//for stabilization
@@ -116,13 +122,18 @@ class scalar_histogram
     std::valarray<double> counter_gathered(0.0, num_data);
     std::valarray<double> sumval_gathered(0.0, num_data);
 
+#ifdef ALPS_HAVE_MPI
     alps_comm.barrier();
     assert(sumval.size()==num_data);
     assert(counter_gathered.size()==num_data);
     assert(sumval_gathered.size()==num_data);
     all_reduce(alps_comm, counter, counter_gathered, std::plus<double>());
     all_reduce(alps_comm, sumval, sumval_gathered, std::plus<double>());
     alps_comm.barrier();
+#else
+    counter_gathered = counter;
+    sumval_gathered = sumval;
+#endif
 
     double maxdist_new = maxdist;
 
@@ -240,11 +251,19 @@ class scalar_histogram_flavors
     return sumval_flavors;
   }
 
+#ifdef ALPS_HAVE_MPI
   double update_cutoff(double cutoff_ratio, double maxdist, double mag, const alps::mpi::communicator& alps_comm) const {
+#else
+  double update_cutoff(double cutoff_ratio, double maxdist, double mag) const {
+#endif
     double maxdist_new = -1.0;
     //for (auto& elem : histograms) {
     for (int ielm=0; ielm<histograms.size(); ++ielm) {
+#ifdef ALPS_HAVE_MPI
       boost::tuple<bool,double> r = histograms[ielm].update_cutoff(cutoff_ratio, maxdist, mag, alps_comm);
+#else
+      boost::tuple<bool,double> r = histograms[ielm].update_cutoff(cutoff_ratio, maxdist, mag);
+#endif
       maxdist_new = std::max(maxdist_new, boost::get<1>(r));
     }
     assert(maxdist_new>0);