Update Lorenz ensemble sample to new OpenMP algebra

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

@@ -18,12 +18,15 @@
 #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;
@@ -33,46 +36,57 @@ struct sys_func {
     const vector<double> &R;
     sys_func( vector<double> &R ) : R(R) {}
 
-    void operator()( const inner_state_type &x , inner_state_type &dxdt , double t , size_t offset ) const {
-        for(size_t i = 0 ; i < x.size() ; i++) {
-            const point_type &xi = x[i];
-            point_type &dxdti = dxdt[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];
+    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() {
-    using namespace boost::numeric::odeint;
+int main(int argc, char **argv) {
+    size_t n = 1024;
+    if(argc > 1) n = boost::lexical_cast<size_t>(argv[1]);
 
-    typedef openmp_state<inner_state_type> state_type;
-
-    const size_t n = 1024;
     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;
 
-    inner_state_type inner(n, point_type(10, 10, 10));
-    state_type state(inner);
+    vector<point_type> inner(n, point_type(10, 10, 10));
+    state_type state;
+    boost::numeric::odeint::copy(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(),
-        openmp_wrapper( sys_func(R) ),
+        sys_func(R),
         state,
         0.0, t_max, dt
     );
 
-    for(size_t i = 0 ; i < state.size() ; i++)
-        copy( state[i].begin(), state[i].end(), ostream_iterator<point_type>(cout, "\n") );
+    boost::numeric::odeint::copy(state, inner);
+    std::copy( inner.begin(), inner.end(), ostream_iterator<point_type>(cout, "\n") );
 
     return 0;
 }

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

@@ -16,7 +16,7 @@
 #include <vector>
 #include <iostream>
 #include <boost/numeric/odeint.hpp>
-#include <boost/numeric/odeint/external/openmp/openmp_algebra.hpp>
+#include <boost/numeric/odeint/external/openmp/openmp.hpp>
 #include "point_type.hpp"
 
 using namespace std;
@@ -61,7 +61,7 @@ int main() {
     typedef runge_kutta4<
         state_type, double,
         state_type, double,
-        basic_openmp_algebra
+        openmp_range_algebra
     > stepper;
 
     const double t_max = 10.0, dt = 0.01;