ceph/src/boost/libs/numeric/odeint/examples/openmp/phase_chain_omp_state.cpp

   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 {
  34     phase_chain_omp_state( double gamma = 0.5 )
  35     : m_gamma( gamma ) { }
  36
  37     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
  62     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
  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 }