diff --git a/CMakeLists.txt b/CMakeLists.txt index 6618e89..35abac9 100755 --- a/CMakeLists.txt +++ b/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) diff --git a/src/impurity.hpp b/src/impurity.hpp index 9f886a0..8ca453d 100755 --- a/src/impurity.hpp +++ b/src/impurity.hpp @@ -39,7 +39,9 @@ #include #include #include -#include +#ifdef ALPS_HAVE_MPI + #include +#endif #include #include @@ -108,7 +110,9 @@ class HybridizationSimulation : public alps::mcbase boost::scoped_ptr 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 diff --git a/src/impurity.ipp b/src/impurity.ipp index e428d73..f034219 100755 --- a/src/impurity.ipp +++ b/src/impurity.ipp @@ -48,7 +48,9 @@ HybridizationSimulation::HybridizationSimulation(parameters_type cons N(static_cast(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::update_MC_parameters() { 1.0/std::max(static_cast(sweeps)/static_cast(interval_update_cutoff),1.0) ); +#ifdef ALPS_HAVE_MPI boost::tuple r_pair = weight_vs_distance.update_cutoff(acc_rate_cutoff, max_distance_pair, mag, comm); +#else + boost::tuple 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 r_shift = weight_vs_distance_shift.update_cutoff(acc_rate_cutoff, max_distance_shift, mag, comm); +#else + boost::tuple 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::max(1,static_cast(BETA/(2.0*max_distance)))); @@ -544,8 +552,7 @@ void HybridizationSimulation::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; } @@ -553,7 +560,7 @@ void HybridizationSimulation::prepare_for_measurement() { 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::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::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_swap0) { 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; } } diff --git a/src/impurity_init.ipp b/src/impurity_init.ipp index 7cca6b8..377a705 100755 --- a/src/impurity_init.ipp +++ b/src/impurity_init.ipp @@ -97,7 +97,7 @@ void HybridizationSimulation::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::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::read_eq_time_two_particle_greens_meas() template void HybridizationSimulation::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; diff --git a/src/main.hpp b/src/main.hpp index a13c1e5..43c4bac 100644 --- a/src/main.hpp +++ b/src/main.hpp @@ -16,20 +16,27 @@ #include "util.hpp" #include -#include #include #include -#include #include - -#include "mympiadapter.hpp" +#ifdef ALPS_HAVE_MPI + #include + #include "mc/mympiadapter.hpp" +#endif + #include "mc/mymcadapter.hpp" #include "postprocess.hpp" #undef BUILD_PYTHON_MODULE +int global_mpi_rank; + template int run_simulation(int argc, const char* argv[], typename alps::parameters_type::type& parameters) { +#ifdef ALPS_HAVE_MPI typedef mymcmpiadapter sim_type; +#else + typedef mymcadapter sim_type; +#endif SOLVER_TYPE::define_parameters(parameters); if (parameters.help_requested(std::cout)) { @@ -37,30 +44,43 @@ int run_simulation(int argc, const char* argv[], typename alps::parameters_type< } char** argv_tmp = const_cast(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 cb = alps::stop_callback(c,size_t(parameters["timelimit"])); +#else + global_mpi_rank = 0; + sim_type sim(parameters); + const boost::function 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::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(results, parameters, ar); +#ifdef ALPS_HAVE_MPI } else{ alps::collect_results(sim); } +#endif + +#ifdef ALPS_HAVE_MPI c.barrier(); +#endif return 0; } diff --git a/src/mc/mymcadapter.hpp b/src/mc/mymcadapter.hpp new file mode 100644 index 0000000..9be6a4d --- /dev/null +++ b/src/mc/mymcadapter.hpp @@ -0,0 +1,57 @@ +#pragma once + +#include +#include +#include + +template 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 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((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(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; +}; diff --git a/src/mympiadapter.hpp b/src/mc/mympiadapter.hpp similarity index 100% rename from src/mympiadapter.hpp rename to src/mc/mympiadapter.hpp diff --git a/src/update_histogram.hpp b/src/update_histogram.hpp index 8e7d786..7bcd412 100755 --- a/src/update_histogram.hpp +++ b/src/update_histogram.hpp @@ -2,6 +2,7 @@ #include +#ifdef ALPS_HAVE_MPI #include /// performs mpi_allreduce() for a valarray of type T. @@ -12,6 +13,7 @@ void all_reduce(const alps::mpi::communicator& comm, const std::valarray& 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(), alps::mpi::is_mpi_op::op(), comm); } +#endif template void rebin(std::valarray& 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 update_cutoff(double cutoff_ratio, double maxdist, double mag, const alps::mpi::communicator& alps_comm) const { +#else + boost::tuple 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,6 +122,7 @@ class scalar_histogram std::valarray counter_gathered(0.0, num_data); std::valarray sumval_gathered(0.0, num_data); +#ifdef ALPS_HAVE_MPI alps_comm.barrier(); assert(sumval.size()==num_data); assert(counter_gathered.size()==num_data); @@ -123,6 +130,10 @@ class scalar_histogram all_reduce(alps_comm, counter, counter_gathered, std::plus()); all_reduce(alps_comm, sumval, sumval_gathered, std::plus()); 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 r = histograms[ielm].update_cutoff(cutoff_ratio, maxdist, mag, alps_comm); +#else + boost::tuple r = histograms[ielm].update_cutoff(cutoff_ratio, maxdist, mag); +#endif maxdist_new = std::max(maxdist_new, boost::get<1>(r)); } assert(maxdist_new>0);