Add an MPI parallelized VEGAS with example

.gitignore

@@ -38,9 +38,10 @@ doc/html
 *.lo
 *.o
 
-# ignore files in m4/ subdir, but track autoconf-archive macros
+# ignore files in m4/ subdir, but track custom macros
 m4/*
 !m4/ax_cxx_compile_stdcxx_11.m4
+!m4/lx_find_mpi.m4
 
 # distribution packages
 *.tar.*

configure.ac

@@ -60,6 +60,10 @@ AC_ARG_ENABLE([examples], AS_HELP_STRING([--enable-examples],
 
 AM_CONDITIONAL([HAVE_EXAMPLES], [test "x$enable_examples" = "xyes"])
 
+AS_IF([test "x$enable_examples" = "xyes"], [
+	LX_FIND_MPI
+])
+
 # add possibility to generate API documentation with Doxygen
 AC_ARG_ENABLE([doxygen], AS_HELP_STRING([--enable-doxygen],
 	[Enable generation of Doxygen documentation]))

examples/Makefile.am

@@ -3,5 +3,8 @@ AM_CPPFLAGS = -I$(top_srcdir)/include
 AM_DEFAULT_SOURCE_EXT = .cpp
 
 noinst_PROGRAMS = \
+	mpi_vegas \
 	simple_vegas \
 	vegas_example
+
+mpi_vegas_LDADD = $(MPI_CXXLDFLAGS)

examples/mpi_vegas.cpp

@@ -0,0 +1,71 @@
+#include <hep/mc.hpp>
+#include <hep/mc-mpi.hpp>
+
+#include <cmath>
+#include <cstddef>
+#include <iostream>
+
+double gauss(hep::mc_point<double> const& sample)
+{
+	double result = 1.0;
+
+	// acos(-1.0) = pi - there is no pi constant in C++!!
+	double factor = 2.0 / std::sqrt(std::acos(-1.0));
+
+	// multiply a (half-)gaussian for every dimension
+	for (std::size_t i = 0; i != sample.point.size(); ++i)
+	{
+		double argument = sample.point[i];
+		result *= std::exp(-argument * argument);
+		result *= factor;
+	}
+
+	return result;
+}
+
+int main(int argc, char* argv[])
+{
+	// initialize MPI
+	MPI_Init(&argc, &argv);
+
+	// get id for this process
+	int rank = 0;
+	MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
+	// dimension of the gauss
+	std::size_t const dimension = 10;
+
+	// five iterations with 10^7 calls each
+	std::vector<std::size_t> iterations(5, 10000000);
+
+	// perform the integration
+	auto results = hep::mpi_vegas<double>(
+		MPI_COMM_WORLD,
+		dimension,
+		iterations,
+		gauss
+	);
+
+	// print numbers in scientific format
+	std::cout.setf(std::ios_base::scientific);
+
+	// print result for iterations
+	if (rank == 0)
+	{
+		// erf is the error function - the integral of the gaussian
+		double reference_result = std::pow(std::erf(1.0), double(dimension));
+
+		std::cout << "reference result is " << reference_result << "\n";
+		std::cout << "results of each iteration:\n";
+		for (std::size_t i = 0; i != results.size(); ++i)
+		{
+			std::cout << i << " (N=" << results[i].calls << ") : I=";
+			std::cout << results[i].value << " +- " << results[i].error << "\n";
+		}
+	}
+
+	// clean up
+	MPI_Finalize();
+
+	return 0;
+}

include/Makefile.am

@@ -2,9 +2,11 @@ nobase_include_HEADERS = \
 	hep/mc/mc_helper.hpp \
 	hep/mc/mc_point.hpp \
 	hep/mc/mc_result.hpp \
+	hep/mc/mpi_vegas.hpp \
 	hep/mc/plain.hpp \
 	hep/mc/vegas.hpp \
-	hep/mc.hpp
+	hep/mc.hpp \
+	hep/mc-mpi.hpp
 
 # create a timestamp file so we can depend on it instead of depending on every
 # single header

include/hep/mc-mpi.hpp

@@ -0,0 +1,33 @@
+#ifndef HEP_MC_MPI_HPP
+#define HEP_MC_MPI_HPP
+
+/*
+ * hep-mc - A Template Library for Monte Carlo Integration
+ * Copyright (C) 2013  Christopher Schwan
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <hep/mc/mpi_vegas.hpp>
+
+namespace hep
+{
+
+/**
+ * \example mpi_vegas.cpp
+ */
+
+}
+
+#endif

include/hep/mc/mpi_vegas.hpp

@@ -0,0 +1,111 @@
+#ifndef HEP_MC_MPI_VEGAS_HPP
+#define HEP_MC_MPI_VEGAS_HPP
+
+#include <hep/mc/vegas.hpp>
+
+#include <cstddef>
+#include <random>
+#include <type_traits>
+#include <vector>
+
+#include <mpi.h>
+
+namespace hep
+{
+
+/**
+ * \ingroup algorithms
+ *
+ * Implements the MPI-parallelized VEGAS algorithm. See \ref vegas() for a more
+ * detailed description on the VEGAS algorithm. In contrast to the
+ * single-threaded versions this function makes sure that every random number
+ * generator is seeded differently so every MPI process yields an independent
+ * result.
+ *
+ * \see The file mpi_vegas.cpp provides an example
+ *
+ * \param communicator The MPI communicator that is used to communicate between
+ *        the different MPI processes.
+ * \param dimensions The number of parameters `function` accepts.
+ * \param iteration_calls The number of function calls that are used to obtain
+ *        a result for each iteration. `iteration_calls.size()` determines the
+ *        number of iterations.
+ * \param function The function that will be integrated over the hypercube. See
+ *        \ref integrands for further explanation.
+ * \param bins The number of bins that the grid will contain for each dimension.
+ * \param alpha The \f$ \alpha \f$ parameter of VEGAS. This parameter should be
+ *        between `1` and `2`.
+ * \param generator The random number generator that will be used to generate
+ *        random points from the hypercube. This generator is not seeded.
+ */
+template <typename T, typename F, typename R = std::mt19937>
+std::vector<vegas_iteration_result<T>> mpi_vegas(
+	MPI_Comm communicator,
+	std::size_t dimensions,
+	std::vector<std::size_t> const& iteration_calls,
+	F&& function,
+	std::size_t bins = 128,
+	T alpha = T(1.5),
+	R&& generator = std::mt19937()
+) {
+	int rank = 0;
+	MPI_Comm_rank(communicator, &rank);
+	int world = 0;
+	MPI_Comm_size(communicator, &world);
+
+	MPI_Datatype type;
+	if (std::is_same<float, T>::value) {
+		type = MPI_FLOAT;
+	} else if (std::is_same<double, T>::value) {
+		type = MPI_DOUBLE;
+	} else if (std::is_same<long double, T>::value) {
+		type = MPI_LONG_DOUBLE;
+	} else {
+		// TODO: error !?
+	}
+
+	// seed every random number generator differently
+	std::size_t r = rank * 10;
+	std::seed_seq sequence{r+0, r+1, r+2, r+3, r+4, r+5, r+6, r+7, r+8, r+9};
+	generator.seed(sequence);
+
+	// create a fresh grid
+	auto grid = vegas_grid<T>(dimensions, bins);
+
+	// vector holding all iteration results
+	std::vector<vegas_iteration_result<T>> results;
+	results.reserve(iteration_calls.size());
+
+	// perform iterations
+	for (auto i = iteration_calls.begin(); i != iteration_calls.end(); ++i)
+	{
+		std::size_t const calls = (*i / world) + (rank < (*i % world) ? 1 : 0);
+		auto result = vegas_iteration(calls, *i, grid, function, generator);
+
+		// add up results
+		MPI_Allreduce(
+			MPI_IN_PLACE,
+			&(result.adjustment_data[0]),
+			result.adjustment_data.size(),
+			type,
+			MPI_SUM,
+			communicator
+		);
+
+		// extract accumulated sum and sum of squares
+		T const& sum = result.adjustment_data[dimensions * bins];
+		T const& sum_of_squares = result.adjustment_data[dimensions * bins + 1];
+
+		// calculate accumulated results
+		results.push_back(vegas_iteration_result<T>(*i, grid,
+			result.adjustment_data));
+
+		grid = vegas_adjust_grid(alpha, grid, result.adjustment_data);
+	}
+
+	return results;
+}
+
+}
+
+#endif