Trivial algebra for OpenMP #6

boost/numeric/odeint/external/openmp/openmp_algebra.hpp

@@ -0,0 +1,65 @@
+/*
+ [auto_generated]
+ boost/numeric/odeint/external/openmp/openmp_algebra.hpp
+
+ [begin_description]
+ Nested parallelized algebra for OpenMP.
+ [end_description]
+
+ Copyright 2009-2011 Karsten Ahnert
+ Copyright 2009-2011 Mario Mulansky
+
+ 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)
+ */
+
+
+#ifndef BOOST_NUMERIC_ODEINT_EXTERNAL_OPENMP_OPENMP_ALGEBRA_HPP_INCLUDED
+#define BOOST_NUMERIC_ODEINT_EXTERNAL_OPENMP_OPENMP_ALGEBRA_HPP_INCLUDED
+
+#include <boost/numeric/odeint/algebra/norm_result_type.hpp>
+#include <boost/numeric/odeint/util/n_ary_helper.hpp>
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+/** \brief OpenMP-parallelized algebra.
+ *
+ * Requires `s.size()` and `s[n]`, i.e. a Random Access Container.
+ */
+struct openmp_algebra
+{
+
+// FIXME: _Pragma is C++11.
+#define OPENMP_ALGEBRA(n) \
+    const size_t len = s0.size(); \
+    _Pragma("omp parallel for") \
+    for( size_t i = 0 ; i < len ; i++ ) \
+        op( BOOST_PP_ENUM_BINARY_PARAMS(n, s, [i] BOOST_PP_INTERCEPT) );
+BOOST_ODEINT_GEN_FOR_EACH(OPENMP_ALGEBRA)
+#undef OPENMP_ALGEBRA
+
+
+    template< class S >
+    static typename norm_result_type< S >::type norm_inf( const S &s )
+    {
+        using std::max;
+        using std::abs;
+        typedef typename norm_result_type< S >::type result_type;
+        result_type init = static_cast< result_type >( 0 );
+#       pragma omp parallel for reduction(max: init)
+        for( size_t i = 0 ; i < s.size() ; ++i )
+            init = max( init , abs(s[i]) );
+        return init;
+    }
+
+};
+
+
+}
+}
+}
+
+#endif

libs/numeric/odeint/examples/Jamfile.v2

@@ -39,3 +39,4 @@ exe bind_member_functions_cpp11 : bind_member_functions_cpp11.cpp : <cxxflags>-s
 # build-project mtl ;
 # build-project ublas ;
 # build-project gmpxx ;
+build-project openmp ;

libs/numeric/odeint/examples/openmp/Jamfile.v2

@@ -0,0 +1,18 @@
+# Copyright 2009 Karsten Ahnert and Mario Mulansky.
+# 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)
+
+project
+    : requirements
+      <include>../../../../..
+      <include>..
+      <define>BOOST_ALL_NO_LIB=1
+      <cxxflags>--std=c++11
+      <toolset>gcc:<cxxflags>-fopenmp
+      <toolset>gcc:<linkflags>-fopenmp
+      <toolset>intel:<cxxflags>-openmp
+      <toolset>intel:<linkflags>-openmp
+    ;
+
+exe lorenz_ensemble : lorenz_ensemble.cpp ;

libs/numeric/odeint/examples/openmp/lorenz_ensemble.cpp

@@ -0,0 +1,96 @@
+/* Boost libs/numeric/odeint/examples/openmp/lorenz_ensemble.cpp
+
+ Copyright 2009-2012 Karsten Ahnert
+ Copyright 2009-2012 Mario Mulansky
+
+ 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)
+ */
+
+#if _OMP
+
+#include <omp.h>
+#include <vector>
+#include <iostream>
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/external/openmp/openmp_algebra.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));
+
+    typedef runge_kutta4<
+        state_type, double,
+        state_type, double,
+        openmp_algebra
+    > stepper;
+
+    const double t_max = 10.0, dt = 0.01;
+
+    const int thr = omp_get_max_threads();
+    for(size_t i = 0 ; i < 10 ; i++) {
+        omp_set_num_threads(thr);
+        const double start0 = omp_get_wtime();
+        integrate_const(
+            stepper(),
+            sys_func(R), X,
+            0.0, t_max, dt);
+        const double delta0 = omp_get_wtime() - start0;
+
+        omp_set_num_threads(1);
+        const double start1 = omp_get_wtime();
+        integrate_const(
+            stepper(),
+            sys_func(R), X,
+            0.0, t_max, dt);
+        const double delta1 = omp_get_wtime() - start1;
+
+        cout << thr << "t=" << delta0 << "s\t1t=" << delta1 << "s\tspeedup=" << (delta1 / delta0) << endl;
+    }
+    return 0;
+}
+
+#else
+
+int main() { return -1; }
+
+#endif