From 2277080c9f63bc4551dfe66fd2dbd5941fbc61b8 Mon Sep 17 00:00:00 2001 From: Pascal Germroth Date: Sat, 15 Jun 2013 18:07:16 +0200 Subject: [PATCH] Trivial algebra for OpenMP #6 --- .../odeint/external/openmp/openmp_algebra.hpp | 65 +++++++++++++ libs/numeric/odeint/examples/Jamfile.v2 | 1 + .../numeric/odeint/examples/openmp/Jamfile.v2 | 18 ++++ .../examples/openmp/lorenz_ensemble.cpp | 96 +++++++++++++++++++ 4 files changed, 180 insertions(+) create mode 100644 boost/numeric/odeint/external/openmp/openmp_algebra.hpp create mode 100644 libs/numeric/odeint/examples/openmp/Jamfile.v2 create mode 100644 libs/numeric/odeint/examples/openmp/lorenz_ensemble.cpp diff --git a/boost/numeric/odeint/external/openmp/openmp_algebra.hpp b/boost/numeric/odeint/external/openmp/openmp_algebra.hpp new file mode 100644 index 00000000..333c18be --- /dev/null +++ b/boost/numeric/odeint/external/openmp/openmp_algebra.hpp @@ -0,0 +1,65 @@ +/* + [auto_generated] + boost/numeric/odeint/external/openmp/openmp_algebra.hpp + + [begin_description] + Nested parallelized algebra for OpenMP. + [end_description] + + Copyright 2009-2011 Karsten Ahnert + Copyright 2009-2011 Mario Mulansky + + Distributed under the Boost Software License, Version 1.0. + (See accompanying file LICENSE_1_0.txt or + copy at http://www.boost.org/LICENSE_1_0.txt) + */ + + +#ifndef BOOST_NUMERIC_ODEINT_EXTERNAL_OPENMP_OPENMP_ALGEBRA_HPP_INCLUDED +#define BOOST_NUMERIC_ODEINT_EXTERNAL_OPENMP_OPENMP_ALGEBRA_HPP_INCLUDED + +#include +#include + +namespace boost { +namespace numeric { +namespace odeint { + +/** \brief OpenMP-parallelized algebra. + * + * Requires `s.size()` and `s[n]`, i.e. a Random Access Container. + */ +struct openmp_algebra +{ + +// FIXME: _Pragma is C++11. +#define OPENMP_ALGEBRA(n) \ + const size_t len = s0.size(); \ + _Pragma("omp parallel for") \ + for( size_t i = 0 ; i < len ; i++ ) \ + op( BOOST_PP_ENUM_BINARY_PARAMS(n, s, [i] BOOST_PP_INTERCEPT) ); +BOOST_ODEINT_GEN_FOR_EACH(OPENMP_ALGEBRA) +#undef OPENMP_ALGEBRA + + + template< class S > + static typename norm_result_type< S >::type norm_inf( const S &s ) + { + using std::max; + using std::abs; + typedef typename norm_result_type< S >::type result_type; + result_type init = static_cast< result_type >( 0 ); +# pragma omp parallel for reduction(max: init) + for( size_t i = 0 ; i < s.size() ; ++i ) + init = max( init , abs(s[i]) ); + return init; + } + +}; + + +} +} +} + +#endif diff --git a/libs/numeric/odeint/examples/Jamfile.v2 b/libs/numeric/odeint/examples/Jamfile.v2 index bad21735..7b9b7ef7 100644 --- a/libs/numeric/odeint/examples/Jamfile.v2 +++ b/libs/numeric/odeint/examples/Jamfile.v2 @@ -39,3 +39,4 @@ exe bind_member_functions_cpp11 : bind_member_functions_cpp11.cpp : -s # build-project mtl ; # build-project ublas ; # build-project gmpxx ; +build-project openmp ; diff --git a/libs/numeric/odeint/examples/openmp/Jamfile.v2 b/libs/numeric/odeint/examples/openmp/Jamfile.v2 new file mode 100644 index 00000000..16c38c2a --- /dev/null +++ b/libs/numeric/odeint/examples/openmp/Jamfile.v2 @@ -0,0 +1,18 @@ +# Copyright 2009 Karsten Ahnert and Mario Mulansky. +# Distributed under the Boost Software License, Version 1.0. (See +# accompanying file LICENSE_1_0.txt or copy at +# http://www.boost.org/LICENSE_1_0.txt) + +project + : requirements + ../../../../.. + .. + BOOST_ALL_NO_LIB=1 + --std=c++11 + gcc:-fopenmp + gcc:-fopenmp + intel:-openmp + intel:-openmp + ; + +exe lorenz_ensemble : lorenz_ensemble.cpp ; diff --git a/libs/numeric/odeint/examples/openmp/lorenz_ensemble.cpp b/libs/numeric/odeint/examples/openmp/lorenz_ensemble.cpp new file mode 100644 index 00000000..04017d02 --- /dev/null +++ b/libs/numeric/odeint/examples/openmp/lorenz_ensemble.cpp @@ -0,0 +1,96 @@ +/* Boost libs/numeric/odeint/examples/openmp/lorenz_ensemble.cpp + + Copyright 2009-2012 Karsten Ahnert + Copyright 2009-2012 Mario Mulansky + + Parallelized Lorenz ensembles + + Distributed under the Boost Software License, Version 1.0. +(See accompanying file LICENSE_1_0.txt or + copy at http://www.boost.org/LICENSE_1_0.txt) + */ + +#if _OMP + +#include +#include +#include +#include +#include +#include "point_type.hpp" + +using namespace std; + +typedef vector vector_type; +typedef point point_type; +typedef vector state_type; + +const double sigma = 10.0; +const double b = 8.0 / 3.0; + +struct sys_func { + const vector_type &R; + sys_func( const vector_type &_R ) : R( _R ) { } + + void operator()( const state_type &x , state_type &dxdt , double t ) const { + const size_t n = x.size(); +# pragma omp parallel for + for(size_t i = 0 ; i < n ; i++) { + const point_type &xi = x[i]; + point_type &dxdti = dxdt[i]; + dxdti[0] = -sigma * (xi[0] - xi[1]); + dxdti[1] = R[i] * xi[0] - xi[1] - xi[0] * xi[2]; + dxdti[2] = -b * xi[2] + xi[0] * xi[1]; + } + } +}; + + +int main() { + using namespace boost::numeric::odeint; + + const size_t n = 1024; + vector_type R(n); + const double Rmin = 0.1, Rmax = 50.0; +# pragma omp parallel for + for(size_t i = 0 ; i < n ; i++) + R[i] = Rmin + (Rmax - Rmin) / (n - 1) * i; + + state_type X(n, point_type(10, 10, 10)); + + typedef runge_kutta4< + state_type, double, + state_type, double, + openmp_algebra + > stepper; + + const double t_max = 10.0, dt = 0.01; + + const int thr = omp_get_max_threads(); + for(size_t i = 0 ; i < 10 ; i++) { + omp_set_num_threads(thr); + const double start0 = omp_get_wtime(); + integrate_const( + stepper(), + sys_func(R), X, + 0.0, t_max, dt); + const double delta0 = omp_get_wtime() - start0; + + omp_set_num_threads(1); + const double start1 = omp_get_wtime(); + integrate_const( + stepper(), + sys_func(R), X, + 0.0, t_max, dt); + const double delta1 = omp_get_wtime() - start1; + + cout << thr << "t=" << delta0 << "s\t1t=" << delta1 << "s\tspeedup=" << (delta1 / delta0) << endl; + } + return 0; +} + +#else + +int main() { return -1; } + +#endif