initial odeint commit

[SVN r81491]

performance/const_step_iterator.cpp

@@ -0,0 +1,164 @@
+/*
+ * odeint_rk4_lorenz.cpp
+ *
+ * Copyright 2009-2012 Karsten Ahnert
+ * Copyright 2009-2012 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)
+ */
+
+#include <algorithm>
+#include <array>
+
+#include <boost/range/algorithm.hpp>
+#include <boost/range/adaptor/filtered.hpp>
+#include <boost/range/numeric.hpp>
+#include <boost/timer.hpp>
+
+#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
+#include <boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp>
+#include <boost/numeric/odeint/stepper/generation.hpp>
+#include <boost/numeric/odeint/iterator/const_step_iterator.hpp>
+#include <boost/numeric/odeint/iterator/const_step_time_iterator.hpp>
+#include <boost/numeric/odeint/integrate/integrate_const.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+typedef std::array< double , 3 > state_type;
+typedef std::vector< double  > state_type2;
+
+std::ostream& operator<<( std::ostream &out , const state_type &x )
+{
+    out << x[0] << "\t" << x[1] << "\t" << x[2];
+    return out;
+}
+
+std::ostream& operator<<( std::ostream &out , const state_type2 &x )
+{
+    if( !x.empty() ) out << x[0];
+    for( size_t i=1 ; i<x.size() ; ++i ) out << "\t" << x[i];
+    return out;
+}
+
+
+const double sigma = 10.0;
+const double R = 28.0;
+const double b = 8.0 / 3.0;
+
+struct lorenz
+{
+    template< class State , class Deriv >
+    void operator()( const State &x , Deriv &dxdt , double t ) const
+    {
+        dxdt[0] = sigma * ( x[1] - x[0] );
+        dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+        dxdt[2] = -b * x[2] + x[0] * x[1];
+    }
+};
+
+boost::timer timer;
+
+void test1_iterator( double t_end = 1000000.0 , double dt = 0.01 )
+{
+    cout << "Test 1 Iterator with array : ";
+    timer.restart();
+
+    runge_kutta4< state_type > stepper;
+    state_type x = {{ 10.0 , 10.0 , 10.0 }};
+    auto first = make_const_step_iterator_begin( stepper , lorenz() , x , 0.0 , t_end , dt );
+    auto last = make_const_step_iterator_end( stepper , lorenz() , x );
+    for( ; first != last ; )
+        ++first;
+
+    cout << "\t" << x << " elapsed time : " << timer.elapsed() << " s" << endl;
+}
+
+void test1_time_iterator( double t_end = 1000000.0 , double dt = 0.01 )
+{
+    cout << "Test 1 Time iterator with array : ";
+    timer.restart();
+
+    runge_kutta4< state_type > stepper;
+    state_type x = {{ 10.0 , 10.0 , 10.0 }};
+    auto first = make_const_step_time_iterator_begin( stepper , lorenz() , x , 0.0 , t_end , dt );
+    auto last = make_const_step_time_iterator_end( stepper , lorenz() , x );
+    for( ; first != last ; )
+        ++first;
+
+    cout << "\t" << x << " elapsed time : " << timer.elapsed() << " s" << endl;
+}
+
+void test1_integrate( double t_end = 1000000.0 , double dt = 0.01 )
+{
+    cout << "Test 1 Integrate with array : ";
+    timer.restart();
+
+    runge_kutta4< state_type > stepper;
+    state_type x = {{ 10.0 , 10.0 , 10.0 }};
+    integrate_const( stepper , lorenz() , x , 0.0 , t_end , dt );
+
+    cout << "\t" << x << " elapsed time : " << timer.elapsed() << " s" << endl;
+}
+
+
+
+
+void test2_iterator( double t_end = 1000000.0 , double dt = 0.01 )
+{
+    cout << "Test 2 Iterator with vector : ";
+    timer.restart();
+
+    runge_kutta4< state_type2 > stepper;
+    state_type2 x = { 10.0 , 10.0 , 10.0 };
+    auto first = make_const_step_iterator_begin( stepper , lorenz() , x , 0.0 , t_end , dt );
+    auto last = make_const_step_iterator_end( stepper , lorenz() , x );
+    for( ; first != last ; )
+        ++first;
+
+    cout << "\t" << x << " elapsed time : " << timer.elapsed() << " s" << endl;
+}
+
+void test2_time_iterator( double t_end = 1000000.0 , double dt = 0.01 )
+{
+    cout << "Test 2 Time iterator with vector : ";
+    timer.restart();
+
+    runge_kutta4< state_type2 > stepper;
+    state_type2 x = { 10.0 , 10.0 , 10.0 };
+    auto first = make_const_step_time_iterator_begin( stepper , lorenz() , x , 0.0 , t_end , dt );
+    auto last = make_const_step_time_iterator_end( stepper , lorenz() , x );
+    for( ; first != last ; )
+        ++first;
+
+    cout << "\t" << x << " elapsed time : " << timer.elapsed() << " s" << endl;
+}
+
+void test2_integrate( double t_end = 1000000.0 , double dt = 0.01 )
+{
+    cout << "Test 2 Integrate with vector : ";
+    timer.restart();
+    runge_kutta4< state_type2 > stepper;
+    state_type2 x = { 10.0 , 10.0 , 10.0 };
+    integrate_const( stepper , lorenz() , x , 0.0 , t_end , dt );
+    cout << "\t" << x << " elapsed time : " << timer.elapsed() << " s" << endl;
+}
+
+
+
+int main( int argc , char **argv )
+{
+    test1_iterator( 1000000.0 );
+    test1_time_iterator( 1000000.0 );
+    test1_integrate( 1000000.0 );
+
+    test2_iterator( 1000000.0 );
+    test2_time_iterator( 1000000.0 );
+    test2_integrate( 1000000.0 );
+
+
+    
+    return 0;
+}