diff options
Diffstat (limited to 'src/boost/libs/numeric/odeint/examples/openmp')
7 files changed, 513 insertions, 0 deletions
diff --git a/src/boost/libs/numeric/odeint/examples/openmp/Jamfile.v2 b/src/boost/libs/numeric/odeint/examples/openmp/Jamfile.v2 new file mode 100644 index 000000000..cef4f670c --- /dev/null +++ b/src/boost/libs/numeric/odeint/examples/openmp/Jamfile.v2 @@ -0,0 +1,23 @@ +# Copyright 2011-2013 Mario Mulansky +# Copyright 2012 Karsten Ahnert +# Copyright 2013 Pascal Germroth +# 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) + +use-project /boost : $(BOOST_ROOT) ; +import openmp : * ; + +project + : requirements + <include>.. + <define>BOOST_ALL_NO_LIB=1 + <library>/boost//timer + [ openmp ] + ; + +exe lorenz_ensemble : lorenz_ensemble.cpp ; +exe lorenz_ensemble_simple : lorenz_ensemble_simple.cpp ; +exe lorenz_ensemble_nested : lorenz_ensemble_nested.cpp ; +exe phase_chain : phase_chain.cpp ; +exe phase_chain_omp_state : phase_chain_omp_state.cpp ; diff --git a/src/boost/libs/numeric/odeint/examples/openmp/lorenz_ensemble.cpp b/src/boost/libs/numeric/odeint/examples/openmp/lorenz_ensemble.cpp new file mode 100644 index 000000000..6717a505b --- /dev/null +++ b/src/boost/libs/numeric/odeint/examples/openmp/lorenz_ensemble.cpp @@ -0,0 +1,91 @@ +/* Boost libs/numeric/odeint/examples/openmp/lorenz_ensemble.cpp + + Copyright 2013 Karsten Ahnert + Copyright 2013 Mario Mulansky + Copyright 2013 Pascal Germroth + + 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) + */ + +#include <omp.h> +#include <vector> +#include <iostream> +#include <iterator> +#include <boost/numeric/odeint.hpp> +#include <boost/numeric/odeint/external/openmp/openmp.hpp> +#include <boost/lexical_cast.hpp> +#include "point_type.hpp" + +using namespace std; +using namespace boost::numeric::odeint; + +typedef point<double, 3> point_type; +typedef vector< point_type > inner_state_type; +typedef openmp_state<point_type> state_type; + +const double sigma = 10.0; +const double b = 8.0 / 3.0; + + +struct sys_func { + const vector<double> &R; + sys_func( vector<double> &R ) : R(R) {} + + void operator()( const state_type &x , state_type &dxdt , double t ) const { +# pragma omp parallel for + for(size_t j = 0 ; j < x.size() ; j++) { + size_t offset = 0; + for(size_t i = 0 ; i < j ; i++) + offset += x[i].size(); + + for(size_t i = 0 ; i < x[j].size() ; i++) { + const point_type &xi = x[j][i]; + point_type &dxdti = dxdt[j][i]; + dxdti[0] = -sigma * (xi[0] - xi[1]); + dxdti[1] = R[offset + i] * xi[0] - xi[1] - xi[0] * xi[2]; + dxdti[2] = -b * xi[2] + xi[0] * xi[1]; + } + } + } +}; + + +int main(int argc, char **argv) { + size_t n = 1024; + if(argc > 1) n = boost::lexical_cast<size_t>(argv[1]); + + vector<double> 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; + + vector<point_type> inner(n, point_type(10, 10, 10)); + state_type state; + split(inner, state); + + cerr << "openmp_state split " << n << " into"; + for(size_t i = 0 ; i != state.size() ; i++) + cerr << ' ' << state[i].size(); + cerr << endl; + + typedef runge_kutta4< state_type, double > stepper; + + const double t_max = 10.0, dt = 0.01; + + integrate_const( + stepper(), + sys_func(R), + state, + 0.0, t_max, dt + ); + + unsplit(state, inner); + std::copy( inner.begin(), inner.end(), ostream_iterator<point_type>(cout, "\n") ); + + return 0; +} diff --git a/src/boost/libs/numeric/odeint/examples/openmp/lorenz_ensemble_nested.cpp b/src/boost/libs/numeric/odeint/examples/openmp/lorenz_ensemble_nested.cpp new file mode 100644 index 000000000..4609c47a4 --- /dev/null +++ b/src/boost/libs/numeric/odeint/examples/openmp/lorenz_ensemble_nested.cpp @@ -0,0 +1,75 @@ +/* Boost libs/numeric/odeint/examples/openmp/lorenz_ensemble_nested.cpp + + Copyright 2013 Karsten Ahnert + Copyright 2013 Pascal Germroth + Copyright 2013 Mario Mulansky + + Parallelized Lorenz ensembles using nested omp algebra + + 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) + */ + +#include <omp.h> +#include <vector> +#include <iostream> +#include <iterator> +#include <boost/numeric/odeint.hpp> +#include <boost/numeric/odeint/external/openmp/openmp.hpp> +#include <boost/lexical_cast.hpp> +#include "point_type.hpp" + +using namespace std; +using namespace boost::numeric::odeint; + +typedef point<double, 3> point_type; +typedef vector< point_type > state_type; + +const double sigma = 10.0; +const double b = 8.0 / 3.0; + + +struct sys_func { + const vector<double> &R; + sys_func( vector<double> &R ) : R(R) {} + + void operator()( const state_type &x , state_type &dxdt , double t ) const { +# pragma omp parallel for + for(size_t i = 0 ; i < x.size() ; i++) { + dxdt[i][0] = -sigma * (x[i][0] - x[i][1]); + dxdt[i][1] = R[i] * x[i][0] - x[i][1] - x[i][0] * x[i][2]; + dxdt[i][2] = -b * x[i][2] + x[i][0] * x[i][1]; + } + } +}; + + +int main(int argc, char **argv) { + size_t n = 1024; + if(argc > 1) n = boost::lexical_cast<size_t>(argv[1]); + + vector<double> 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 state( n , point_type(10, 10, 10) ); + + typedef runge_kutta4< state_type, double , state_type , double , + openmp_nested_algebra<vector_space_algebra> > stepper; + + const double t_max = 10.0, dt = 0.01; + + integrate_const( + stepper(), + sys_func(R), + state, + 0.0, t_max, dt + ); + + std::copy( state.begin(), state.end(), ostream_iterator<point_type>(cout, "\n") ); + + return 0; +} diff --git a/src/boost/libs/numeric/odeint/examples/openmp/lorenz_ensemble_simple.cpp b/src/boost/libs/numeric/odeint/examples/openmp/lorenz_ensemble_simple.cpp new file mode 100644 index 000000000..fbaf24999 --- /dev/null +++ b/src/boost/libs/numeric/odeint/examples/openmp/lorenz_ensemble_simple.cpp @@ -0,0 +1,79 @@ +/* Boost libs/numeric/odeint/examples/openmp/lorenz_ensemble_simple.cpp + + Copyright 2013 Karsten Ahnert + Copyright 2013 Mario Mulansky + Copyright 2013 Pascal Germroth + + 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) + */ + +#include <omp.h> +#include <vector> +#include <iostream> +#include <boost/numeric/odeint.hpp> +#include <boost/numeric/odeint/external/openmp/openmp.hpp> +#include "point_type.hpp" + +using namespace std; + +typedef vector<double> vector_type; +typedef point<double, 3> point_type; +typedef vector<point_type> 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)); + + const double t_max = 10.0, dt = 0.01; + + // Simple stepper with constant step size + // typedef runge_kutta4<state_type, double, state_type, double, + // openmp_range_algebra> stepper; + + // integrate_const(stepper(), sys_func(R), X, 0.0, t_max, dt); + + // Controlled stepper with variable step size + typedef runge_kutta_fehlberg78<state_type, double, state_type, double, + openmp_range_algebra> error_stepper_type; + typedef controlled_runge_kutta<error_stepper_type> controlled_stepper_type; + controlled_stepper_type controlled_stepper; + + integrate_adaptive(controlled_stepper, sys_func(R), X, 0.0, t_max, dt); + + copy( X.begin(), X.end(), ostream_iterator<point_type>(cout, "\n") ); + + return 0; +} diff --git a/src/boost/libs/numeric/odeint/examples/openmp/openmp.jam b/src/boost/libs/numeric/odeint/examples/openmp/openmp.jam new file mode 100644 index 000000000..226bcf577 --- /dev/null +++ b/src/boost/libs/numeric/odeint/examples/openmp/openmp.jam @@ -0,0 +1,52 @@ +# Copyright 2013 Karsten Ahnert +# Copyright 2013 Mario Mulansky +# Copyright 2013 Pascal Germroth +# 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) + + +# Only builds target with supported OpenMP enabled toolsets. +# +# use as: +# exe omp : omp.cpp : [ openmp ] ; +# +rule openmp return + # default + <build>no + # GNU C++ + <toolset>gcc:<cxxflags>-fopenmp + <toolset>gcc:<linkflags>-fopenmp + <toolset>gcc:<build>yes + # Microsoft Visual C++ + <toolset>msvc:<cxxflags>/openmp + <toolset>msvc:<linkflags>/openmp + <toolset>msvc:<build>yes + # Intel C++ + <toolset>intel-linux:<cxxflags>-openmp + <toolset>intel-linux:<linkflags>-openmp + <toolset>intel-linux:<build>yes + <toolset>intel-win:<cxxflags>-Qopenmp + <toolset>intel-win:<linkflags>-Qopenmp + <toolset>intel-win:<build>yes + # HP aC++ + <toolset>acc:<cxxflags>+Oopenmp + <toolset>acc:<linkflags>+Oopenmp + <toolset>acc:<build>yes + # Sun Studio + <toolset>sun:<cxxflags>-xopenmp + <toolset>sun:<linkflags>-xopenmp + <toolset>sun:<build>yes + # IBM XL + <toolset>vacpp:<cxxflags>-qsmp=omp + <toolset>vacpp:<linkflags>-qsmp=omp + <toolset>vacpp:<build>yes + # PG++ + <toolset>pgi:<cxxflags>-mp + <toolset>pgi:<linkflags>-mp + <toolset>pgi:<build>yes + # Pathscale + <toolset>pathscale:<cxxflags>-mp + <toolset>pathscale:<linkflags>-mp + <toolset>pathscale:<build>yes + ; diff --git a/src/boost/libs/numeric/odeint/examples/openmp/phase_chain.cpp b/src/boost/libs/numeric/odeint/examples/openmp/phase_chain.cpp new file mode 100644 index 000000000..08876fcce --- /dev/null +++ b/src/boost/libs/numeric/odeint/examples/openmp/phase_chain.cpp @@ -0,0 +1,95 @@ +/* + * phase_chain.cpp + * + * Example of OMP parallelization with odeint + * + * Copyright 2013 Karsten Ahnert + * Copyright 2013 Mario Mulansky + * Copyright 2013 Pascal Germroth + * 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) + */ + +#include <iostream> +#include <vector> +#include <boost/random.hpp> +#include <boost/timer/timer.hpp> +//[phase_chain_openmp_header +#include <omp.h> +#include <boost/numeric/odeint.hpp> +#include <boost/numeric/odeint/external/openmp/openmp.hpp> +//] + +using namespace std; +using namespace boost::numeric::odeint; +using boost::timer::cpu_timer; +using boost::math::double_constants::pi; + +//[phase_chain_vector_state +typedef std::vector< double > state_type; +//] + +//[phase_chain_rhs +struct phase_chain +{ + phase_chain( double gamma = 0.5 ) + : m_gamma( gamma ) { } + + void operator()( const state_type &x , state_type &dxdt , double /* t */ ) const + { + const size_t N = x.size(); + #pragma omp parallel for schedule(runtime) + for(size_t i = 1 ; i < N - 1 ; ++i) + { + dxdt[i] = coupling_func( x[i+1] - x[i] ) + + coupling_func( x[i-1] - x[i] ); + } + dxdt[0 ] = coupling_func( x[1 ] - x[0 ] ); + dxdt[N-1] = coupling_func( x[N-2] - x[N-1] ); + } + + double coupling_func( double x ) const + { + return sin( x ) - m_gamma * ( 1.0 - cos( x ) ); + } + + double m_gamma; +}; +//] + + +int main( int argc , char **argv ) +{ + //[phase_chain_init + size_t N = 131101; + state_type x( N ); + boost::random::uniform_real_distribution<double> distribution( 0.0 , 2.0*pi ); + boost::random::mt19937 engine( 0 ); + generate( x.begin() , x.end() , boost::bind( distribution , engine ) ); + //] + + //[phase_chain_stepper + typedef runge_kutta4< + state_type , double , + state_type , double , + openmp_range_algebra + > stepper_type; + //] + + //[phase_chain_scheduling + int chunk_size = N/omp_get_max_threads(); + omp_set_schedule( omp_sched_static , chunk_size ); + //] + + cpu_timer timer; + //[phase_chain_integrate + integrate_n_steps( stepper_type() , phase_chain( 1.2 ) , + x , 0.0 , 0.01 , 100 ); + //] + double run_time = static_cast<double>(timer.elapsed().wall) * 1.0e-9; + std::cerr << run_time << "s" << std::endl; + // copy(x.begin(), x.end(), ostream_iterator<double>(cout, "\n")); + + return 0; +} diff --git a/src/boost/libs/numeric/odeint/examples/openmp/phase_chain_omp_state.cpp b/src/boost/libs/numeric/odeint/examples/openmp/phase_chain_omp_state.cpp new file mode 100644 index 000000000..8627f58ad --- /dev/null +++ b/src/boost/libs/numeric/odeint/examples/openmp/phase_chain_omp_state.cpp @@ -0,0 +1,98 @@ +/* + * phase_chain_omp_state.cpp + * + * Example of OMP parallelization with odeint + * + * Copyright 2013 Karsten Ahnert + * Copyright 2013 Mario Mulansky + * Copyright 2013 Pascal Germroth + * 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) + */ + +#include <iostream> +#include <vector> +#include <boost/random.hpp> +#include <boost/timer/timer.hpp> +#include <omp.h> +#include <boost/numeric/odeint.hpp> +#include <boost/numeric/odeint/external/openmp/openmp.hpp> + +#include <boost/numeric/odeint.hpp> + +using namespace std; +using namespace boost::numeric::odeint; +using boost::timer::cpu_timer; +using boost::math::double_constants::pi; + +typedef openmp_state<double> state_type; + +//[phase_chain_state_rhs +struct phase_chain_omp_state +{ + phase_chain_omp_state( double gamma = 0.5 ) + : m_gamma( gamma ) { } + + void operator()( const state_type &x , state_type &dxdt , double /* t */ ) const + { + const size_t N = x.size(); + #pragma omp parallel for schedule(runtime) + for(size_t n = 0 ; n < N ; ++n) + { + const size_t M = x[n].size(); + for(size_t m = 1 ; m < M-1 ; ++m) + { + dxdt[n][m] = coupling_func( x[n][m+1] - x[n][m] ) + + coupling_func( x[n][m-1] - x[n][m] ); + } + dxdt[n][0] = coupling_func( x[n][1] - x[n][0] ); + if( n > 0 ) + { + dxdt[n][0] += coupling_func( x[n-1].back() - x[n].front() ); + } + dxdt[n][M-1] = coupling_func( x[n][M-2] - x[n][M-1] ); + if( n < N-1 ) + { + dxdt[n][M-1] += coupling_func( x[n+1].front() - x[n].back() ); + } + } + } + + double coupling_func( double x ) const + { + return sin( x ) - m_gamma * ( 1.0 - cos( x ) ); + } + + double m_gamma; +}; +//] + + +int main( int argc , char **argv ) +{ + //[phase_chain_state_init + const size_t N = 131101; + vector<double> x( N ); + boost::random::uniform_real_distribution<double> distribution( 0.0 , 2.0*pi ); + boost::random::mt19937 engine( 0 ); + generate( x.begin() , x.end() , boost::bind( distribution , engine ) ); + const size_t blocks = omp_get_max_threads(); + state_type x_split( blocks ); + split( x , x_split ); + //] + + + cpu_timer timer; + //[phase_chain_state_integrate + integrate_n_steps( runge_kutta4<state_type>() , phase_chain_omp_state( 1.2 ) , + x_split , 0.0 , 0.01 , 100 ); + unsplit( x_split , x ); + //] + + double run_time = static_cast<double>(timer.elapsed().wall) * 1.0e-9; + std::cerr << run_time << "s" << std::endl; + // copy(x.begin(), x.end(), ostream_iterator<double>(cout, "\n")); + + return 0; +} |