From 82df16841bbec41bd4acfacea900eb8e54d485c8 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Tue, 11 Jun 2013 17:16:21 +0200 Subject: [PATCH] Add an MPI parallelized VEGAS with example --- .gitignore | 3 +- configure.ac | 4 + examples/Makefile.am | 3 + examples/mpi_vegas.cpp | 71 ++++++++++++ include/Makefile.am | 4 +- include/hep/mc-mpi.hpp | 33 ++++++ include/hep/mc/mpi_vegas.hpp | 111 +++++++++++++++++++ m4/lx_find_mpi.m4 | 204 +++++++++++++++++++++++++++++++++++ 8 files changed, 431 insertions(+), 2 deletions(-) create mode 100644 examples/mpi_vegas.cpp create mode 100644 include/hep/mc-mpi.hpp create mode 100644 include/hep/mc/mpi_vegas.hpp create mode 100644 m4/lx_find_mpi.m4 diff --git a/.gitignore b/.gitignore index c9f9a3a..9535bfc 100644 --- a/.gitignore +++ b/.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.* diff --git a/configure.ac b/configure.ac index d837268..fee9358 100644 --- a/configure.ac +++ b/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])) diff --git a/examples/Makefile.am b/examples/Makefile.am index a8ecfbf..8d6c8d9 100644 --- a/examples/Makefile.am +++ b/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) diff --git a/examples/mpi_vegas.cpp b/examples/mpi_vegas.cpp new file mode 100644 index 0000000..c3e55bb --- /dev/null +++ b/examples/mpi_vegas.cpp @@ -0,0 +1,71 @@ +#include +#include + +#include +#include +#include + +double gauss(hep::mc_point 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 iterations(5, 10000000); + + // perform the integration + auto results = hep::mpi_vegas( + 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; +} diff --git a/include/Makefile.am b/include/Makefile.am index 3e3b1b9..0c9458c 100644 --- a/include/Makefile.am +++ b/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 diff --git a/include/hep/mc-mpi.hpp b/include/hep/mc-mpi.hpp new file mode 100644 index 0000000..b473241 --- /dev/null +++ b/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 . + */ + +#include + +namespace hep +{ + +/** + * \example mpi_vegas.cpp + */ + +} + +#endif diff --git a/include/hep/mc/mpi_vegas.hpp b/include/hep/mc/mpi_vegas.hpp new file mode 100644 index 0000000..e6e7f3c --- /dev/null +++ b/include/hep/mc/mpi_vegas.hpp @@ -0,0 +1,111 @@ +#ifndef HEP_MC_MPI_VEGAS_HPP +#define HEP_MC_MPI_VEGAS_HPP + +#include + +#include +#include +#include +#include + +#include + +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 +std::vector> mpi_vegas( + MPI_Comm communicator, + std::size_t dimensions, + std::vector 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::value) { + type = MPI_FLOAT; + } else if (std::is_same::value) { + type = MPI_DOUBLE; + } else if (std::is_same::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(dimensions, bins); + + // vector holding all iteration results + std::vector> 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(*i, grid, + result.adjustment_data)); + + grid = vegas_adjust_grid(alpha, grid, result.adjustment_data); + } + + return results; +} + +} + +#endif diff --git a/m4/lx_find_mpi.m4 b/m4/lx_find_mpi.m4 new file mode 100644 index 0000000..7e9ac9b --- /dev/null +++ b/m4/lx_find_mpi.m4 @@ -0,0 +1,204 @@ +################################################################################################# +# Copyright (c) 2010, Lawrence Livermore National Security, LLC. +# Produced at the Lawrence Livermore National Laboratory +# Written by Todd Gamblin, tgamblin@llnl.gov. +# LLNL-CODE-417602 +# All rights reserved. +# +# This file is part of Libra. For details, see http://github.com/tgamblin/libra. +# Please also read the LICENSE file for further information. +# +# Redistribution and use in source and binary forms, with or without modification, are +# permitted provided that the following conditions are met: +# +# * Redistributions of source code must retain the above copyright notice, this list of +# conditions and the disclaimer below. +# * Redistributions in binary form must reproduce the above copyright notice, this list of +# conditions and the disclaimer (as noted below) in the documentation and/or other materials +# provided with the distribution. +# * Neither the name of the LLNS/LLNL nor the names of its contributors may be used to endorse +# or promote products derived from this software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS +# OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +# MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL +# LAWRENCE LIVERMORE NATIONAL SECURITY, LLC, THE U.S. DEPARTMENT OF ENERGY OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, +# WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +################################################################################################# + +# +# LX_FIND_MPI() +# ------------------------------------------------------------------------ +# This macro finds an MPI compiler and extracts includes and libraries from +# it for use in automake projects. The script exports the following variables: +# +# AC_DEFINE variables: +# HAVE_MPI AC_DEFINE'd to 1 if we found MPI +# +# AC_SUBST variables: +# MPICC Name of MPI compiler +# MPI_CFLAGS Includes and defines for MPI C compilation +# MPI_CLDFLAGS Libraries and library paths for linking MPI C programs +# +# MPICXX Name of MPI C++ compiler +# MPI_CXXFLAGS Includes and defines for MPI C++ compilation +# MPI_CXXLDFLAGS Libraries and library paths for linking MPI C++ programs +# +# MPIF77 Name of MPI Fortran 77 compiler +# MPI_F77FLAGS Includes and defines for MPI Fortran 77 compilation +# MPI_F77LDFLAGS Libraries and library paths for linking MPI Fortran 77 programs +# +# MPIFC Name of MPI Fortran compiler +# MPI_FFLAGS Includes and defines for MPI Fortran compilation +# MPI_FLDFLAGS Libraries and library paths for linking MPI Fortran programs +# +# Shell variables output by this macro: +# have_C_mpi 'yes' if we found MPI for C, 'no' otherwise +# have_CXX_mpi 'yes' if we found MPI for C++, 'no' otherwise +# have_F77_mpi 'yes' if we found MPI for F77, 'no' otherwise +# have_F_mpi 'yes' if we found MPI for Fortran, 'no' otherwise +# +AC_DEFUN([LX_FIND_MPI], +[ + AC_LANG_CASE( + [C], [ + AC_REQUIRE([AC_PROG_CC]) + if [[ ! -z "$MPICC" ]]; then + LX_QUERY_MPI_COMPILER(MPICC, [$MPICC], C) + else + LX_QUERY_MPI_COMPILER(MPICC, [mpicc mpiicc mpixlc mpipgcc], C) + fi + ], + [C++], [ + AC_REQUIRE([AC_PROG_CXX]) + if [[ ! -z "$MPICXX" ]]; then + LX_QUERY_MPI_COMPILER(MPICXX, [$MPICXX], CXX) + else + LX_QUERY_MPI_COMPILER(MPICXX, [mpicxx mpiCC mpic++ mpig++ mpiicpc mpipgCC mpixlC], CXX) + fi + ], + [F77], [ + AC_REQUIRE([AC_PROG_F77]) + if [[ ! -z "$MPIF77" ]]; then + LX_QUERY_MPI_COMPILER(MPIF77, [$MPIF77], F77) + else + LX_QUERY_MPI_COMPILER(MPIF77, [mpif77 mpiifort mpixlf77 mpixlf77_r], F77) + fi + ], + [Fortran], [ + AC_REQUIRE([AC_PROG_FC]) + if [[ ! -z "$MPIFC" ]]; then + LX_QUERY_MPI_COMPILER(MPIFC, [$MPIFC], F) + else + mpi_default_fc="mpif95 mpif90 mpigfortran mpif2003" + mpi_intel_fc="mpiifort" + mpi_xl_fc="mpixlf95 mpixlf95_r mpixlf90 mpixlf90_r mpixlf2003 mpixlf2003_r" + mpi_pg_fc="mpipgf95 mpipgf90" + LX_QUERY_MPI_COMPILER(MPIFC, [$mpi_default_fc $mpi_intel_fc $mpi_xl_fc $mpi_pg_fc], F) + fi + ]) +]) + + +# +# LX_QUERY_MPI_COMPILER([compiler-var-name], [compiler-names], [output-var-prefix]) +# ------------------------------------------------------------------------ +# AC_SUBST variables: +# MPI_FLAGS Includes and defines for MPI compilation +# MPI_LDFLAGS Libraries and library paths for linking MPI C programs +# +# Shell variables output by this macro: +# found_mpi_flags 'yes' if we were able to get flags, 'no' otherwise +# +AC_DEFUN([LX_QUERY_MPI_COMPILER], +[ + # Try to find a working MPI compiler from the supplied names + AC_PATH_PROGS($1, [$2], [not-found]) + + # Figure out what the compiler responds to to get it to show us the compile + # and link lines. After this part of the macro, we'll have a valid + # lx_mpi_command_line + echo -n "Checking whether $$1 responds to '-showme:compile'... " + lx_mpi_compile_line=`$$1 -showme:compile 2>/dev/null` + if [[ "$?" -eq 0 ]]; then + echo yes + lx_mpi_link_line=`$$1 -showme:link 2>/dev/null` + else + echo no + echo -n "Checking whether $$1 responds to '-showme'... " + lx_mpi_command_line=`$$1 -showme 2>/dev/null` + if [[ "$?" -ne 0 ]]; then + echo no + echo -n "Checking whether $$1 responds to '-compile-info'... " + lx_mpi_compile_line=`$$1 -compile-info 2>/dev/null` + if [[ "$?" -eq 0 ]]; then + echo yes + lx_mpi_link_line=`$$1 -link-info 2>/dev/null` + else + echo no + echo -n "Checking whether $$1 responds to '-show'... " + lx_mpi_command_line=`$$1 -show 2>/dev/null` + if [[ "$?" -eq 0 ]]; then + echo yes + else + echo no + fi + fi + else + echo yes + fi + fi + + if [[ ! -z "$lx_mpi_compile_line" -a ! -z "$lx_mpi_link_line" ]]; then + lx_mpi_command_line="$lx_mpi_compile_line $lx_mpi_link_line" + fi + + if [[ ! -z "$lx_mpi_command_line" ]]; then + # Now extract the different parts of the MPI command line. Do these separately in case we need to + # parse them all out in future versions of this macro. + lx_mpi_defines=` echo "$lx_mpi_command_line" | grep -o -- '\(^\| \)-D\([[^\"[:space:]]]\+\|\"[[^\"[:space:]]]\+\"\)'` + lx_mpi_includes=` echo "$lx_mpi_command_line" | grep -o -- '\(^\| \)-I\([[^\"[:space:]]]\+\|\"[[^\"[:space:]]]\+\"\)'` + lx_mpi_link_paths=` echo "$lx_mpi_command_line" | grep -o -- '\(^\| \)-L\([[^\"[:space:]]]\+\|\"[[^\"[:space:]]]\+\"\)'` + lx_mpi_libs=` echo "$lx_mpi_command_line" | grep -o -- '\(^\| \)-l\([[^\"[:space:]]]\+\|\"[[^\"[:space:]]]\+\"\)'` + lx_mpi_link_args=` echo "$lx_mpi_command_line" | grep -o -- '\(^\| \)-Wl,\([[^\"[:space:]]]\+\|\"[[^\"[:space:]]]\+\"\)'` + + # Create variables and clean up newlines and multiple spaces + MPI_$3FLAGS="$lx_mpi_defines $lx_mpi_includes" + MPI_$3LDFLAGS="$lx_mpi_link_paths $lx_mpi_libs $lx_mpi_link_args" + MPI_$3FLAGS=` echo "$MPI_$3FLAGS" | tr '\n' ' ' | sed 's/^[[ \t]]*//;s/[[ \t]]*$//' | sed 's/ +/ /g'` + MPI_$3LDFLAGS=`echo "$MPI_$3LDFLAGS" | tr '\n' ' ' | sed 's/^[[ \t]]*//;s/[[ \t]]*$//' | sed 's/ +/ /g'` + + OLD_CPPFLAGS=$CPPFLAGS + OLD_LIBS=$LIBS + CPPFLAGS=$MPI_$3FLAGS + LIBS=$MPI_$3LDFLAGS + + AC_TRY_LINK([#include ], + [int rank, size; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &size);], + [# Add a define for testing at compile time. + AC_DEFINE([HAVE_MPI], [1], [Define to 1 if you have MPI libs and headers.]) + have_$3_mpi='yes'], + [# zero out mpi flags so we don't link against the faulty library. + MPI_$3FLAGS="" + MPI_$3LDFLAGS="" + have_$3_mpi='no']) + + # AC_SUBST everything. + AC_SUBST($1) + AC_SUBST(MPI_$3FLAGS) + AC_SUBST(MPI_$3LDFLAGS) + + LIBS=$OLD_LIBS + CPPFLAGS=$OLD_CPPFLAGS + else + Echo Unable to find suitable MPI Compiler. Try setting $1. + have_$3_mpi='no' + fi +]) +