Squashed 'third_party/boostorg/odeint/' content from commit 6ff2719
Change-Id: If4892e29c1a5e6cf3a7aa51486a2725c251b0c7d
git-subtree-dir: third_party/boostorg/odeint
git-subtree-split: 6ff2719b6907b86596c3d43e88c1bcfdf29df560
diff --git a/examples/openmp/Jamfile.v2 b/examples/openmp/Jamfile.v2
new file mode 100644
index 0000000..cef4f67
--- /dev/null
+++ b/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/examples/openmp/lorenz_ensemble.cpp b/examples/openmp/lorenz_ensemble.cpp
new file mode 100644
index 0000000..6717a50
--- /dev/null
+++ b/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/examples/openmp/lorenz_ensemble_nested.cpp b/examples/openmp/lorenz_ensemble_nested.cpp
new file mode 100644
index 0000000..4609c47
--- /dev/null
+++ b/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/examples/openmp/lorenz_ensemble_simple.cpp b/examples/openmp/lorenz_ensemble_simple.cpp
new file mode 100644
index 0000000..fbaf249
--- /dev/null
+++ b/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/examples/openmp/openmp.jam b/examples/openmp/openmp.jam
new file mode 100644
index 0000000..226bcf5
--- /dev/null
+++ b/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/examples/openmp/phase_chain.cpp b/examples/openmp/phase_chain.cpp
new file mode 100644
index 0000000..08876fc
--- /dev/null
+++ b/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/examples/openmp/phase_chain_omp_state.cpp b/examples/openmp/phase_chain_omp_state.cpp
new file mode 100644
index 0000000..8627f58
--- /dev/null
+++ b/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;
+}