1 /*
2 * phase_chain_omp_state.cpp
3 *
4 * Example of OMP parallelization with odeint
5 *
6 * Copyright 2013 Karsten Ahnert
7 * Copyright 2013 Mario Mulansky
8 * Copyright 2013 Pascal Germroth
9 * Distributed under the Boost Software License, Version 1.0. (See
10 * accompanying file LICENSE_1_0.txt or copy at
11 * http://www.boost.org/LICENSE_1_0.txt)
12 */
13
14 #include <iostream>
15 #include <vector>
16 #include <boost/random.hpp>
17 #include <boost/timer/timer.hpp>
18 #include <omp.h>
19 #include <boost/numeric/odeint.hpp>
20 #include <boost/numeric/odeint/external/openmp/openmp.hpp>
21
22 #include <boost/numeric/odeint.hpp>
23
24 using namespace std;
25 using namespace boost::numeric::odeint;
26 using boost::timer::cpu_timer;
27 using boost::math::double_constants::pi;
28
29 typedef openmp_state<double> state_type;
30
31 //[phase_chain_state_rhs
32 struct phase_chain_omp_state
33 {
phase_chain_omp_statephase_chain_omp_state34 phase_chain_omp_state( double gamma = 0.5 )
35 : m_gamma( gamma ) { }
36
operator ()phase_chain_omp_state37 void operator()( const state_type &x , state_type &dxdt , double /* t */ ) const
38 {
39 const size_t N = x.size();
40 #pragma omp parallel for schedule(runtime)
41 for(size_t n = 0 ; n < N ; ++n)
42 {
43 const size_t M = x[n].size();
44 for(size_t m = 1 ; m < M-1 ; ++m)
45 {
46 dxdt[n][m] = coupling_func( x[n][m+1] - x[n][m] ) +
47 coupling_func( x[n][m-1] - x[n][m] );
48 }
49 dxdt[n][0] = coupling_func( x[n][1] - x[n][0] );
50 if( n > 0 )
51 {
52 dxdt[n][0] += coupling_func( x[n-1].back() - x[n].front() );
53 }
54 dxdt[n][M-1] = coupling_func( x[n][M-2] - x[n][M-1] );
55 if( n < N-1 )
56 {
57 dxdt[n][M-1] += coupling_func( x[n+1].front() - x[n].back() );
58 }
59 }
60 }
61
coupling_funcphase_chain_omp_state62 double coupling_func( double x ) const
63 {
64 return sin( x ) - m_gamma * ( 1.0 - cos( x ) );
65 }
66
67 double m_gamma;
68 };
69 //]
70
71
main(int argc,char ** argv)72 int main( int argc , char **argv )
73 {
74 //[phase_chain_state_init
75 const size_t N = 131101;
76 vector<double> x( N );
77 boost::random::uniform_real_distribution<double> distribution( 0.0 , 2.0*pi );
78 boost::random::mt19937 engine( 0 );
79 generate( x.begin() , x.end() , boost::bind( distribution , engine ) );
80 const size_t blocks = omp_get_max_threads();
81 state_type x_split( blocks );
82 split( x , x_split );
83 //]
84
85
86 cpu_timer timer;
87 //[phase_chain_state_integrate
88 integrate_n_steps( runge_kutta4<state_type>() , phase_chain_omp_state( 1.2 ) ,
89 x_split , 0.0 , 0.01 , 100 );
90 unsplit( x_split , x );
91 //]
92
93 double run_time = static_cast<double>(timer.elapsed().wall) * 1.0e-9;
94 std::cerr << run_time << "s" << std::endl;
95 // copy(x.begin(), x.end(), ostream_iterator<double>(cout, "\n"));
96
97 return 0;
98 }
99