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

   1 /*
   2  * phase_chain.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 //[phase_chain_openmp_header
  19 #include <omp.h>
  20 #include <boost/numeric/odeint.hpp>
  21 #include <boost/numeric/odeint/external/openmp/openmp.hpp>
  22 //]
  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 //[phase_chain_vector_state
  30 typedef std::vector< double > state_type;
  31 //]
  32
  33 //[phase_chain_rhs
  34 struct phase_chain
  35 {
  36     phase_chain( double gamma = 0.5 )
  37     : m_gamma( gamma ) { }
  38
  39     void operator()( const state_type &x , state_type &dxdt , double /* t */ ) const
  40     {
  41         const size_t N = x.size();
  42         #pragma omp parallel for schedule(runtime)
  43         for(size_t i = 1 ; i < N - 1 ; ++i)
  44         {
  45             dxdt[i] = coupling_func( x[i+1] - x[i] ) +
  46                       coupling_func( x[i-1] - x[i] );
  47         }
  48         dxdt[0  ] = coupling_func( x[1  ] - x[0  ] );
  49         dxdt[N-1] = coupling_func( x[N-2] - x[N-1] );
  50     }
  51
  52     double coupling_func( double x ) const
  53     {
  54         return sin( x ) - m_gamma * ( 1.0 - cos( x ) );
  55     }
  56
  57     double m_gamma;
  58 };
  59 //]
  60
  61
  62 int main( int argc , char **argv )
  63 {
  64     //[phase_chain_init
  65     size_t N = 131101;
  66     state_type x( N );
  67     boost::random::uniform_real_distribution<double> distribution( 0.0 , 2.0*pi );
  68     boost::random::mt19937 engine( 0 );
  69     generate( x.begin() , x.end() , boost::bind( distribution , engine ) );
  70     //]
  71
  72     //[phase_chain_stepper
  73     typedef runge_kutta4<
  74                       state_type , double ,
  75                       state_type , double ,
  76                       openmp_range_algebra
  77                     > stepper_type;
  78     //]
  79
  80     //[phase_chain_scheduling
  81     int chunk_size = N/omp_get_max_threads();
  82     omp_set_schedule( omp_sched_static , chunk_size );
  83     //]
  84
  85     cpu_timer timer;
  86     //[phase_chain_integrate
  87     integrate_n_steps( stepper_type() , phase_chain( 1.2 ) ,
  88                        x , 0.0 , 0.01 , 100 );
  89     //]
  90     double run_time = static_cast<double>(timer.elapsed().wall) * 1.0e-9;
  91     std::cerr << run_time << "s" << std::endl;
  92     // copy(x.begin(), x.end(), ostream_iterator<double>(cout, "\n"));
  93
  94     return 0;
  95 }