Subject: [Boost-commit] svn:boost r71871 - in sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta: . performance
From: mario.mulansky_at_[hidden]
Date: 2011-05-11 11:10:47
Author: mariomulansky
Date: 2011-05-11 11:10:44 EDT (Wed, 11 May 2011)
New Revision: 71871
URL: http://svn.boost.org/trac/boost/changeset/71871
Log:
reorganized test cases
performance tests are killing me!!! it seems we get higher performance when using structs instead of functions for the rhs...
Added:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk4_lorenz.cpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk4_lorenz.cpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk4_lorenz.cpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_lorenz.cpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_lorenz_def_alg.cpp
- copied, changed from r71869, /sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_def_alg.cpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/rk_performance_test_case.hpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/rt_generic_rk4_lorenz.cpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/rt_algebra.hpp
- copied unchanged from r71869, /sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/rt_algebra.hpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/rt_explicit_rk.hpp (contents, props changed)
Removed:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_def_alg.cpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/rt_algebra.hpp
Text files modified:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_explicit_rk_new.hpp | 3
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/Jamfile | 28 ++++++++++-
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk4.cpp | 23 +++++----
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk4.cpp | 1
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk4.cpp | 26 ++++++----
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4.cpp | 9 ++-
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_lorenz_def_alg.cpp | 96 +++++++++++++++++----------------------
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/perf_analysis.py | 8 +-
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/perf_tests.py | 12 +++-
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/results.dat | 2
10 files changed, 116 insertions(+), 92 deletions(-)
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_explicit_rk_new.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_explicit_rk_new.hpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_explicit_rk_new.hpp 2011-05-11 11:10:44 EDT (Wed, 11 May 2011)
@@ -16,6 +16,7 @@
#include <boost/mpl/size_t.hpp>
#include <boost/fusion/container.hpp>
+#include <boost/fusion/algorithm/iteration.hpp>
#include <boost/array.hpp>
@@ -193,7 +194,7 @@
template< class System >
- void inline do_step( System &system , state_type &x , const double t , const double dt )
+ void inline do_step( System system , state_type &x , const double t , const double dt )
{
fusion::for_each( m_stages , calculate_stage< System >( system , x , m_x_tmp , m_F , t , dt ) );
}
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/Jamfile (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/Jamfile 2011-05-11 11:10:44 EDT (Wed, 11 May 2011)
@@ -22,27 +22,49 @@
: <toolset>intel:<cxxflags>-inline-forceinline
;
+exe generic_rk4_lorenz
+ : generic_rk4_lorenz.cpp
+ : <toolset>intel:<cxxflags>-inline-forceinline
+ ;
+
exe odeint_rk4
- : odeint_rk4.cpp
+ : odeint_rk4.cpp
+# : <cxxflags>-Winline
+ ;
+
+exe odeint_rk4_lorenz
+ : odeint_rk4_lorenz.cpp
;
-exe odeint_rk4_def_alg
- : odeint_rk4_def_alg.cpp
+exe odeint_rk4_lorenz_def_alg
+ : odeint_rk4_lorenz_def_alg.cpp
;
exe nr_rk4
: nr_rk4.cpp
;
+
+exe nr_rk4_lorenz
+ : nr_rk4_lorenz.cpp
+ ;
+
exe rt_generic_rk4
: rt_generic_rk4.cpp
;
+
+exe rt_generic_rk4_lorenz
+ : rt_generic_rk4_lorenz.cpp
+ ;
exe gsl_rk4
: gsl_rk4.cpp libgsl libgslcblas
;
+exe gsl_rk4_lorenz
+ : gsl_rk4_lorenz.cpp libgsl libgslcblas
+ ;
exe generic_rk54ck
: generic_rk54ck.cpp
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk4.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk4.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk4.cpp 2011-05-11 11:10:44 EDT (Wed, 11 May 2011)
@@ -39,17 +39,18 @@
typedef boost::array< double , 3 > state_type;
typedef explicit_rk< state_type , 4 > rk4_fusion_type;
-
-inline void lorenz( const state_type &x , state_type &dxdt , const double t )
+struct lorenz
{
- const double sigma = 10.0;
- const double R = 28.0;
- const double b = 8.0 / 3.0;
- dxdt[0] = sigma * ( x[1] - x[0] );
- dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
- dxdt[2] = x[0]*x[1] - b * x[2];
-}
-
+ inline void operator()( const state_type &x , state_type &dxdt , const double t )
+ {
+ const double sigma = 10.0;
+ const double R = 28.0;
+ const double b = 8.0 / 3.0;
+ dxdt[0] = sigma * ( x[1] - x[0] );
+ dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+ dxdt[2] = x[0]*x[1] - b * x[2];
+ }
+};
const size_t loops = 20;
@@ -85,7 +86,7 @@
timer.restart();
for( size_t i=0 ; i<num_of_steps ; ++i, t+=dt )
- rk4_fusion.do_step( lorenz , x , t , dt );
+ rk4_fusion.do_step( lorenz() , x , t , dt );
acc( timer.elapsed() );
clog.precision( 3 );
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk4_lorenz.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk4_lorenz.cpp 2011-05-11 11:10:44 EDT (Wed, 11 May 2011)
@@ -0,0 +1,72 @@
+#include <boost/array.hpp>
+
+#include "../fusion_explicit_rk_new.hpp"
+
+#include "rk_performance_test_case.hpp"
+
+typedef boost::array< double , 3 > state_type;
+typedef explicit_rk< state_type , 4 > rk4_fusion_type;
+
+struct lorenz
+{
+ inline void operator()( const state_type &x , state_type &dxdt , const double t ) const
+ {
+ const double sigma = 10.0;
+ const double R = 28.0;
+ const double b = 8.0 / 3.0;
+ dxdt[0] = sigma * ( x[1] - x[0] );
+ dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+ dxdt[2] = x[0]*x[1] - b * x[2];
+ }
+};
+
+typedef rk4_fusion_type::coef_a_type coef_a_type;
+typedef rk4_fusion_type::coef_b_type coef_b_type;
+typedef rk4_fusion_type::coef_c_type coef_c_type;
+
+const boost::array< double , 1 > a1 = {{ 0.5 }};
+const boost::array< double , 2 > a2 = {{ 0.0 , 0.5 }};
+const boost::array< double , 3 > a3 = {{ 0.0 , 0.0 , 1.0 }};
+
+const coef_a_type a = fusion::make_vector( a1 , a2 , a3 );
+const coef_b_type b = {{ 1.0/6 , 1.0/3 , 1.0/3 , 1.0/6 }};
+const coef_c_type c = {{ 0.0 , 0.5 , 0.5 , 1.0 }};
+
+class fusion_wrapper
+{
+
+public:
+
+ fusion_wrapper() : m_stepper( a , b , c )
+ { }
+
+ void reset_init_cond()
+ {
+ m_x[0] = 10.0 * rand() / RAND_MAX;
+ m_x[1] = 10.0 * rand() / RAND_MAX;
+ m_x[2] = 10.0 * rand() / RAND_MAX;
+ m_t = 0.0;
+ }
+
+ inline void do_step( const double dt )
+ {
+ m_stepper.do_step( lorenz() , m_x , m_t , dt );
+ }
+
+ double state( const size_t i ) const
+ { return m_x[i]; }
+
+private:
+ state_type m_x;
+ double m_t;
+ rk4_fusion_type m_stepper;
+};
+
+
+
+int main()
+{
+ fusion_wrapper stepper;
+
+ run( stepper );
+}
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk4.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk4.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk4.cpp 2011-05-11 11:10:44 EDT (Wed, 11 May 2011)
@@ -81,6 +81,7 @@
clog << acc << " " << x[0] << endl;
}
cout << acc << endl;
+ gsl_odeiv_step_free( s );
return 0;
}
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk4_lorenz.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk4_lorenz.cpp 2011-05-11 11:10:44 EDT (Wed, 11 May 2011)
@@ -0,0 +1,77 @@
+/*
+ * gsl_rk4_lorenz.cpp
+ *
+ * Created on: May 11, 2011
+ * Author: mario
+ */
+
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_odeiv.h>
+
+#include "rk_performance_test_case.hpp"
+
+const size_t dim = 3;
+
+int lorenz_gsl( const double t , const double y[] , double f[] , void *params)
+{
+ const double sigma = 10.0;
+ const double R = 28.0;
+ const double b = 8.0 / 3.0;
+
+ f[0] = sigma * ( y[1] - y[0] );
+ f[1] = R * y[0] - y[1] - y[0] * y[2];
+ f[2] = y[0]*y[1] - b * y[2];
+ return GSL_SUCCESS;
+}
+
+class gsl_wrapper
+{
+public:
+
+ gsl_wrapper()
+ {
+ m_s = gsl_odeiv_step_alloc( gsl_odeiv_step_rk4 , dim);
+ m_sys.function = lorenz_gsl;
+ m_sys.jacobian = 0;
+ m_sys.dimension = dim;
+ m_sys.params = 0;
+ }
+
+ void reset_init_cond()
+ {
+ m_x[0] = 10.0 * rand() / RAND_MAX;
+ m_x[1] = 10.0 * rand() / RAND_MAX;
+ m_x[2] = 10.0 * rand() / RAND_MAX;
+ m_t = 0.0;
+ }
+
+ inline void do_step( const double dt )
+ {
+ gsl_odeiv_step_apply ( m_s , m_t , dt , m_x , m_x_err , 0 , 0 , &m_sys );
+ //m_t += dt;
+ }
+
+ double state( const size_t i ) const
+ { return m_x[i]; }
+
+ ~gsl_wrapper()
+ {
+ gsl_odeiv_step_free( m_s );
+ }
+
+private:
+ double m_x[dim];
+ double m_x_err[dim];
+ double m_t;
+ gsl_odeiv_step *m_s;
+ gsl_odeiv_system m_sys;
+};
+
+
+
+int main()
+{
+ gsl_wrapper stepper;
+
+ run( stepper , 20000000 / 3 , 1E-10 * 3);
+}
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk4.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk4.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk4.cpp 2011-05-11 11:10:44 EDT (Wed, 11 May 2011)
@@ -36,19 +36,23 @@
typedef boost::array< double , 3 > state_type;
-inline void lorenz( const state_type &x , state_type &dxdt , const double t )
-{
- const double sigma = 10.0;
- const double R = 28.0;
- const double b = 8.0 / 3.0;
- dxdt[0] = sigma * ( x[1] - x[0] );
- dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
- dxdt[2] = x[0]*x[1] - b * x[2];
-}
+struct lorenz {
+
+ static const double sigma = 10.0;
+ static const double R = 28.0;
+ static const double b = 8.0 / 3.0;
+
+ void operator()( const state_type &x , state_type &dxdt , const double t ) const
+ {
+ dxdt[0] = sigma * ( x[1] - x[0] );
+ dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+ dxdt[2] = x[0]*x[1] - b * x[2];
+ }
+};
template< class System , typename T , size_t dim >
-void rk4_step( System &sys , boost::array< T , dim > &x , const double t , const double dt )
+void rk4_step( System sys , boost::array< T , dim > &x , const double t , const double dt )
{ // fast rk4 implementation adapted from the book 'Numerical Recipes'
size_t i;
const double hh = dt*0.5;
@@ -97,7 +101,7 @@
timer.restart();
for( size_t i=0 ; i<num_of_steps ; ++i, t+=dt )
- rk4_step( lorenz , x , t , dt );
+ rk4_step( lorenz() , x , t , dt );
acc( timer.elapsed() );
clog.precision( 3 );
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk4_lorenz.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk4_lorenz.cpp 2011-05-11 11:10:44 EDT (Wed, 11 May 2011)
@@ -0,0 +1,89 @@
+/*
+ * nr_rk4_lorenz.cpp
+ *
+ * Created on: May 11, 2011
+ * Author: mario
+ */
+
+#include <boost/array.hpp>
+
+#include "rk_performance_test_case.hpp"
+
+const size_t dim = 3;
+
+typedef boost::array< double , dim > state_type;
+
+struct lorenz
+{
+ inline void operator()( const state_type &x , state_type &dxdt , const double t ) const
+ {
+ const double sigma = 10.0;
+ const double R = 28.0;
+ const double b = 8.0 / 3.0;
+ dxdt[0] = sigma * ( x[1] - x[0] );
+ dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+ dxdt[2] = x[0]*x[1] - b * x[2];
+ }
+};
+
+template< class System , typename T , size_t dim >
+void rk4_step( const System sys , boost::array< T , dim > &x , const double t , const double dt )
+{ // fast rk4 implementation adapted from the book 'Numerical Recipes'
+ size_t i;
+ const double hh = dt*0.5;
+ const double h6 = dt/6.0;
+ const double th = t+hh;
+ boost::array< T , dim > dydx , dym , dyt , yt;
+
+ sys( x , dydx , t );
+
+ for( i=0 ; i<dim ; i++ )
+ yt[i] = x[i] + hh*dydx[i];
+
+ sys( yt , dyt , th );
+ for( i=0 ; i<dim ; i++ )
+ yt[i] = x[i] + hh*dyt[i];
+
+ sys( yt , dym , th );
+ for( i=0 ; i<dim ; i++ ) {
+ yt[i] = x[i] + dt*dym[i];
+ dym[i] += dyt[i];
+ }
+ sys( yt , dyt , t+dt );
+ for( i=0 ; i<dim ; i++ )
+ x[i] += h6*( dydx[i] + dyt[i] + 2.0*dym[i] );
+}
+
+
+class nr_wrapper
+{
+public:
+ void reset_init_cond()
+ {
+ m_x[0] = 10.0 * rand() / RAND_MAX;
+ m_x[1] = 10.0 * rand() / RAND_MAX;
+ m_x[2] = 10.0 * rand() / RAND_MAX;
+ m_t = 0.0;
+ }
+
+ inline void do_step( const double dt )
+ {
+ rk4_step( lorenz() , m_x , m_t , dt );
+ }
+
+ double state( const size_t i ) const
+ { return m_x[i]; }
+
+private:
+ state_type m_x;
+ double m_t;
+};
+
+
+
+int main()
+{
+ nr_wrapper stepper;
+
+ run( stepper );
+}
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4.cpp 2011-05-11 11:10:44 EDT (Wed, 11 May 2011)
@@ -40,11 +40,12 @@
boost::numeric::odeint::array_algebra > rk4_odeint_type;
-inline void lorenz( const state_type &x , state_type &dxdt , const double t )
+const double sigma = 10.0;
+const double R = 28.0;
+const double b = 8.0 / 3.0;
+
+__inline__ void lorenz( const state_type &x , state_type &dxdt , const double t )
{
- const double sigma = 10.0;
- const double R = 28.0;
- const double b = 8.0 / 3.0;
dxdt[0] = sigma * ( x[1] - x[0] );
dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
dxdt[2] = x[0]*x[1] - b * x[2];
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_def_alg.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_def_alg.cpp 2011-05-11 11:10:44 EDT (Wed, 11 May 2011)
+++ (empty file)
@@ -1,84 +0,0 @@
-/*
- * odeint_rk4.cpp
- *
- * Created on: Apr 28, 2011
- * Author: mario
- */
-
-#include <iostream>
-#include <fstream>
-
-#include <boost/array.hpp>
-
-#include <boost/numeric/odeint/stepper/explicit_rk4.hpp>
-#include <boost/numeric/odeint/algebra/array_algebra.hpp>
-#include <boost/accumulators/accumulators.hpp>
-#include <boost/accumulators/statistics.hpp>
-#include <boost/timer.hpp>
-
-#define tab "\t"
-
-using namespace std;
-using namespace boost::accumulators;
-
-typedef accumulator_set<
- double , stats< tag::mean , tag::variance >
- > accumulator_type;
-
-ostream& operator<<( ostream& out , accumulator_type &acc )
-{
- out << boost::accumulators::mean( acc ) << tab;
-// out << boost::accumulators::variance( acc ) << tab;
- return out;
-}
-
-typedef boost::timer timer_type;
-
-
-typedef boost::array< double , 3 > state_type;
-typedef boost::numeric::odeint::explicit_rk4< state_type > rk4_odeint_type;
-
-inline void lorenz( const state_type &x , state_type &dxdt , const double t )
-{
- const double sigma = 10.0;
- const double R = 28.0;
- const double b = 8.0 / 3.0;
- dxdt[0] = sigma * ( x[1] - x[0] );
- dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
- dxdt[2] = x[0]*x[1] - b * x[2];
-}
-
-
-const size_t loops = 20;
-
-int main( int argc , char **argv )
-{
- rk4_odeint_type rk4_odeint;
-
- const size_t num_of_steps = 20000000;
- const double dt = 1E-10;
-
- accumulator_type acc;
- timer_type timer;
-
- srand( 12312354 );
-
- for( size_t n=0 ; n<loops ; ++n )
- {
- state_type x = {{ 10.0 * rand()/RAND_MAX ,
- 10.0 * rand()/RAND_MAX ,
- 10.0 * rand()/RAND_MAX }};
- double t = 0.0;
-
- timer.restart();
- for( size_t i=0 ; i<num_of_steps ; ++i, t+=dt )
- rk4_odeint.do_step( lorenz , x , t , dt );
- acc( timer.elapsed() );
-
- clog.precision( 3 );
- clog.width( 5 );
- clog << acc << " " << x[0] << endl;
- }
- cout << acc << endl;
- return 0;
-}
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_lorenz.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_lorenz.cpp 2011-05-11 11:10:44 EDT (Wed, 11 May 2011)
@@ -0,0 +1,75 @@
+/*
+ * odeint_rk4_lorenz.cpp
+ *
+ * Created on: May 11, 2011
+ * Author: mario
+ */
+
+#include <boost/array.hpp>
+
+#include <boost/numeric/odeint/stepper/explicit_rk4.hpp>
+#include <boost/numeric/odeint/algebra/array_algebra.hpp>
+
+#include "rk_performance_test_case.hpp"
+
+typedef boost::array< double , 3 > state_type;
+typedef boost::numeric::odeint::explicit_rk4< state_type , double , state_type , double ,
+ boost::numeric::odeint::array_algebra > rk4_odeint_type;
+
+struct lorenz
+{
+ inline void operator()( const state_type &x , state_type &dxdt , const double t ) const
+ {
+ const double sigma = 10.0;
+ const double R = 28.0;
+ const double b = 8.0 / 3.0;
+ dxdt[0] = sigma * ( x[1] - x[0] );
+ dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+ dxdt[2] = x[0]*x[1] - b * x[2];
+ }
+};
+
+inline void lorenz_func( const state_type &x , state_type &dxdt , const double t )
+{
+ const double sigma = 10.0;
+ const double R = 28.0;
+ const double b = 8.0 / 3.0;
+ dxdt[0] = sigma * ( x[1] - x[0] );
+ dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+ dxdt[2] = x[0]*x[1] - b * x[2];
+}
+
+class odeint_wrapper
+{
+public:
+ void reset_init_cond()
+ {
+ m_x[0] = 10.0 * rand() / RAND_MAX;
+ m_x[1] = 10.0 * rand() / RAND_MAX;
+ m_x[2] = 10.0 * rand() / RAND_MAX;
+ m_t = 0.0;
+ }
+
+ inline void do_step( const double dt )
+ {
+ m_stepper.do_step( lorenz() , m_x , m_t , dt );
+ //m_t += dt;
+ }
+
+ double state( const size_t i ) const
+ { return m_x[i]; }
+
+private:
+ state_type m_x;
+ double m_t;
+ rk4_odeint_type m_stepper;
+};
+
+
+
+int main()
+{
+ odeint_wrapper stepper;
+
+ run( stepper );
+}
Copied: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_lorenz_def_alg.cpp (from r71869, /sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_def_alg.cpp)
==============================================================================
--- /sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_def_alg.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_lorenz_def_alg.cpp 2011-05-11 11:10:44 EDT (Wed, 11 May 2011)
@@ -1,44 +1,34 @@
/*
- * odeint_rk4.cpp
+ * odeint_rk4_lorenz_def_alg.cpp
*
* Created on: Apr 28, 2011
* Author: mario
*/
-#include <iostream>
-#include <fstream>
-
#include <boost/array.hpp>
#include <boost/numeric/odeint/stepper/explicit_rk4.hpp>
#include <boost/numeric/odeint/algebra/array_algebra.hpp>
-#include <boost/accumulators/accumulators.hpp>
-#include <boost/accumulators/statistics.hpp>
-#include <boost/timer.hpp>
-
-#define tab "\t"
-
-using namespace std;
-using namespace boost::accumulators;
-
-typedef accumulator_set<
- double , stats< tag::mean , tag::variance >
- > accumulator_type;
-
-ostream& operator<<( ostream& out , accumulator_type &acc )
-{
- out << boost::accumulators::mean( acc ) << tab;
-// out << boost::accumulators::variance( acc ) << tab;
- return out;
-}
-
-typedef boost::timer timer_type;
+#include "rk_performance_test_case.hpp"
typedef boost::array< double , 3 > state_type;
typedef boost::numeric::odeint::explicit_rk4< state_type > rk4_odeint_type;
-inline void lorenz( const state_type &x , state_type &dxdt , const double t )
+struct lorenz
+{
+ inline void operator()( const state_type &x , state_type &dxdt , const double t ) const
+ {
+ const double sigma = 10.0;
+ const double R = 28.0;
+ const double b = 8.0 / 3.0;
+ dxdt[0] = sigma * ( x[1] - x[0] );
+ dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+ dxdt[2] = x[0]*x[1] - b * x[2];
+ }
+};
+
+inline void lorenz_func( const state_type &x , state_type &dxdt , const double t )
{
const double sigma = 10.0;
const double R = 28.0;
@@ -48,37 +38,37 @@
dxdt[2] = x[0]*x[1] - b * x[2];
}
+class odeint_wrapper
+{
+public:
+ void reset_init_cond()
+ {
+ m_x[0] = 10.0 * rand() / RAND_MAX;
+ m_x[1] = 10.0 * rand() / RAND_MAX;
+ m_x[2] = 10.0 * rand() / RAND_MAX;
+ m_t = 0.0;
+ }
-const size_t loops = 20;
+ inline void do_step( const double dt )
+ {
+ m_stepper.do_step( lorenz_func , m_x , m_t , dt );
+ //m_t += dt;
+ }
-int main( int argc , char **argv )
-{
- rk4_odeint_type rk4_odeint;
+ double state( const size_t i ) const
+ { return m_x[i]; }
- const size_t num_of_steps = 20000000;
- const double dt = 1E-10;
+private:
+ state_type m_x;
+ double m_t;
+ rk4_odeint_type m_stepper;
+};
- accumulator_type acc;
- timer_type timer;
- srand( 12312354 );
- for( size_t n=0 ; n<loops ; ++n )
- {
- state_type x = {{ 10.0 * rand()/RAND_MAX ,
- 10.0 * rand()/RAND_MAX ,
- 10.0 * rand()/RAND_MAX }};
- double t = 0.0;
-
- timer.restart();
- for( size_t i=0 ; i<num_of_steps ; ++i, t+=dt )
- rk4_odeint.do_step( lorenz , x , t , dt );
- acc( timer.elapsed() );
-
- clog.precision( 3 );
- clog.width( 5 );
- clog << acc << " " << x[0] << endl;
- }
- cout << acc << endl;
- return 0;
+int main()
+{
+ odeint_wrapper stepper;
+
+ run( stepper );
}
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/perf_analysis.py
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/perf_analysis.py (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/perf_analysis.py 2011-05-11 11:10:44 EDT (Wed, 11 May 2011)
@@ -38,14 +38,14 @@
title("Runge-Kutta 4" , fontsize=20)
bar( arange(5) , means_rk4 , bar_width , color='blue' , linewidth=4 , edgecolor='blue' , yerr = error_rk4 , ecolor='red') #, elinewidth=2, ecolor='red' )
xlim( -0.5 , 4.5+bar_width )
-xticks( arange(5)+bar_width/2 , ('odeint' , 'generic' , 'nr' , 'gsl' , 'rt gen' ) )
+xticks( arange(5)+bar_width/2 , ('odeint' , 'generic' , 'NR' , 'GSL' , 'rt gen' ) )
ylabel('Performance in %' , fontsize=20)
figure(2)
title("Runge-Kutta 5(4) Cash-Karp" , fontsize=20)
-bar( arange(4) , means_rk54ck , bar_width , color='blue' , linewidth=4 , edgecolor='blue' , yerr = error_rk54ck , ecolor='red') #, elinewidth=2, ecolor='red' )
+bar( arange(4) , means_rk54ck , bar_width , color='blue' , linewidth=4 , edgecolor='blue' , yerr = error_rk54ck , ecolor='red' )
xlim( -0.5 , 3.5+bar_width )
-xticks( arange(4)+bar_width/2 , ('odeint' , 'generic' , 'nr' , 'gsl' ) )
+xticks( arange(4)+bar_width/2 , ('odeint' , 'generic' , 'NR' , 'GSL' ) )
ylabel('Performance in %' , fontsize=20)
-show()
\ No newline at end of file
+show()
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/perf_tests.py
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/perf_tests.py (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/perf_tests.py 2011-05-11 11:10:44 EDT (Wed, 11 May 2011)
@@ -1,14 +1,17 @@
from os import popen
+from os import system
from os.path import isfile
-bin_path = "bin/gcc-4.4/release/"
-#bin_path = "bin/intel-linux/release/"
+toolset = "intel"
+
+#bin_path = "bin/gcc-4.3/release/"
+bin_path = "bin/intel-linux/release/"
#bin_path = "bin\\msvc-9.0express\\release\\threading-multi\\"
#extension = ".exe"
extension = ""
-bins = [ "odeint_rk4" , "odeint_rk4_def_alg" , "generic_rk4" , "nr_rk4" , "gsl_rk4" , "rt_generic_rk4" ,
- "odeint_rk54ck" , "odeint_rk54ck_def_alg" , "generic_rk54ck" , "nr_rk54ck" , "gsl_rk54ck" ]
+bins = [ "odeint_rk4_lorenz" , "odeint_rk4_lorenz_def_alg" , "generic_rk4_lorenz" , "nr_rk4_lorenz" , "gsl_rk4_lorenz" , "rt_generic_rk4_lorenz" ]
+# "odeint_rk54ck" , "odeint_rk54ck_def_alg" , "generic_rk54ck" , "nr_rk54ck" , "gsl_rk54ck" ]
results = []
@@ -16,6 +19,7 @@
print
for bin in bins:
+ system( "bjam --toolset=" + toolset + " -a " + bin );
if isfile( bin_path + bin + extension):
print "Running" , bin
res = popen( bin_path+bin+extension ).read()
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/results.dat
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/results.dat (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/results.dat 2011-05-11 11:10:44 EDT (Wed, 11 May 2011)
@@ -8,7 +8,7 @@
gcc 4.4.1 | 0.60 | 0.59 | 0.60 | 0.60 | 1.77 | 2.08 | Opteron 2224 @ 3 GHz
gcc 4.4.1 | 0.731 | 0.7315 | 0.729 | 0.7785 | 1.7755 | 2.8295 | Core 2 Duo P8400 @ 2.26GHz
gcc 4.4.1 | 0.7125 | 0.7095 | 0.7085 | 0.7555 | 2.0065 | 2.915 | Core 2 Quad Q6600 @ 2.40GHz
-gcc 4.5.0 | 0.54 | 0.77 | 0.77 | 0.80 | 1.07 | 1.08 | Core i7 870 @ 2.93 GHz
+gcc 4.5.0 | 0.54 | 0.92 | 0.51 | 0.47 | 1.13 | 1.16 | Core i7 870 @ 2.93 GHz
gcc 4.6.0 | 0.54 | 0.54 | 0.65 | 0.47 | 1.06 | 1.08 | Core i7 870 @ 2.93 GHz
icc 12.0.2 | 0.90 | 0.81 | 0.85 | 1.05 | 1.07 | 1.57 | Core i7 870 @ 2.93 GHz
icc 11.1 | 1.11 | 0.98 | 0.98 | 0.63 | 1.28 | 1.70 | Xeon X5650 @ 2.67 GHz
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/rk_performance_test_case.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/rk_performance_test_case.hpp 2011-05-11 11:10:44 EDT (Wed, 11 May 2011)
@@ -0,0 +1,56 @@
+/*
+ * rk_performance_test_case.hpp
+ *
+ * Created on: May 11, 2011
+ * Author: mario
+ */
+
+#include <iostream>
+#include <boost/accumulators/accumulators.hpp>
+#include <boost/accumulators/statistics.hpp>
+#include <boost/timer.hpp>
+
+#define tab "\t"
+
+using namespace std;
+using namespace boost::accumulators;
+
+typedef accumulator_set<
+ double , stats< tag::mean , tag::variance >
+ > accumulator_type;
+
+ostream& operator<<( ostream& out , accumulator_type &acc )
+{
+ out << boost::accumulators::mean( acc ) << tab;
+// out << boost::accumulators::variance( acc ) << tab;
+ return out;
+}
+
+typedef boost::timer timer_type;
+
+
+template< class Stepper >
+void run( Stepper &stepper , const size_t num_of_steps = 20000000 , const double dt = 1E-10 )
+{
+ const size_t loops = 20;
+
+ accumulator_type acc;
+ timer_type timer;
+
+ srand( 12312354 );
+
+ for( size_t n=0 ; n<loops ; ++n )
+ {
+ stepper.reset_init_cond( );
+
+ timer.restart();
+ for( size_t i = 0 ; i < num_of_steps ; ++i )
+ stepper.do_step( dt );
+ acc(timer.elapsed());
+
+ clog.precision(3);
+ clog.width(5);
+ clog << acc << " " << stepper.state(0) << endl;
+ }
+ cout << acc << endl;
+}
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/rt_algebra.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/rt_algebra.hpp 2011-05-11 11:10:44 EDT (Wed, 11 May 2011)
+++ (empty file)
@@ -1,19 +0,0 @@
-#include <vector>
-
-using namespace std;
-
-struct rt_algebra
-{
- template< typename T , size_t dim >
- inline static void foreach( boost::array< T , dim > & x_tmp , const boost::array< T , dim > &x ,
- const vector< double > &a ,
- const boost::array< T , dim > *k_vector , const double dt )
- {
- for( size_t i=0 ; i<dim ; ++i )
- {
- x_tmp[i] = x[i];
- for( size_t j = 0 ; j<a.size() ; ++j )
- x_tmp[i] += a[j]*dt*k_vector[j][i];
- }
- }
-};
\ No newline at end of file
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/rt_generic_rk4_lorenz.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/rt_generic_rk4_lorenz.cpp 2011-05-11 11:10:44 EDT (Wed, 11 May 2011)
@@ -0,0 +1,82 @@
+/*
+ * rt_generic_rk4_lorenz.cpp
+ *
+ * Created on: May 11, 2011
+ * Author: mario
+ */
+
+#include <boost/array.hpp>
+
+#include "../rt_explicit_rk.hpp"
+
+#include "rk_performance_test_case.hpp"
+
+typedef boost::array< double , 3 > state_type;
+
+typedef rt_explicit_rk< state_type > rk_stepper_type;
+
+const size_t stage_count = 4;
+
+inline void lorenz( const state_type &x , state_type &dxdt , const double t )
+{
+ const double sigma = 10.0;
+ const double R = 28.0;
+ const double b = 8.0 / 3.0;
+ dxdt[0] = sigma * ( x[1] - x[0] );
+ dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+ dxdt[2] = x[0]*x[1] - b * x[2];
+}
+
+
+
+class rt_generic_wrapper
+{
+public:
+
+ rt_generic_wrapper()
+ {
+ rk_stepper_type::coeff_a_type a( stage_count-1 );
+ a[0].resize(1); a[0][0] = 0.5;
+ a[1].resize(2); a[1][0] = 0.0; a[1][1] = 0.5;
+ a[2].resize(3); a[2][0] = 0.0; a[2][1] = 0.0; a[2][2] = 1.0;
+
+ rk_stepper_type::coeff_b_type b( stage_count );
+ b[0] = 1.0/6; b[1] = 1.0/3; b[2] = 1.0/3; b[3] = 1.0/6;
+
+ rk_stepper_type::coeff_c_type c( stage_count );
+ c[0] = 0.0; c[1] = 0.5; c[2] = 0.5; c[3] = 1.0;
+
+ m_stepper = rk_stepper_type( stage_count , a , b , c );
+ }
+
+ void reset_init_cond()
+ {
+ m_x[0] = 10.0 * rand() / RAND_MAX;
+ m_x[1] = 10.0 * rand() / RAND_MAX;
+ m_x[2] = 10.0 * rand() / RAND_MAX;
+ m_t = 0.0;
+ }
+
+ inline void do_step( const double dt )
+ {
+ m_stepper.do_step( lorenz , m_x , m_t , dt );
+ //m_t += dt;
+ }
+
+ double state( const size_t i ) const
+ { return m_x[i]; }
+
+private:
+ state_type m_x;
+ double m_t;
+ rk_stepper_type m_stepper;
+};
+
+
+
+int main()
+{
+ rt_generic_wrapper stepper;
+
+ run( stepper );
+}
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/rt_explicit_rk.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/rt_explicit_rk.hpp 2011-05-11 11:10:44 EDT (Wed, 11 May 2011)
@@ -0,0 +1,69 @@
+/*
+ * rt_explicit_rk.hpp
+ *
+ * Created on: May 11, 2011
+ * Author: mario
+ */
+#ifndef RT_EXPLICIT_RK_HPP_
+#define RT_EXPLICIT_RK_HPP_
+
+#include <vector>
+
+#include "rt_algebra.hpp"
+
+using namespace std;
+
+template< class StateType >
+class rt_explicit_rk
+{
+public:
+ typedef StateType state_type;
+ typedef vector< vector< double > > coeff_a_type;
+ typedef vector< double > coeff_b_type;
+ typedef vector< double > coeff_c_type;
+
+ rt_explicit_rk()
+ {
+ m_F = 0;
+ }
+
+ rt_explicit_rk( size_t stage_count , coeff_a_type &a , coeff_b_type &b , coeff_c_type &c )
+ : m_s( stage_count ) , m_a( a ) , m_b( b ) , m_c( c )
+ {
+ m_F = new state_type[ m_s ];
+ }
+
+ ~rt_explicit_rk() { if( m_F ) delete[] m_F; }
+
+ template< class System >
+ void do_step( System &sys , state_type &x , const double t , const double dt )
+ {
+ // first stage separately
+ sys( x , m_F[0] , t + m_c[0]*t );
+ if( m_s == 1 )
+ rt_algebra::foreach( x , x , m_b , m_F , dt );
+ else
+ rt_algebra::foreach( m_x_tmp , x , m_a[0] , m_F , dt );
+
+ for( size_t stage = 2 ; stage <= m_s ; ++stage )
+ {
+ sys( m_x_tmp , m_F[stage-1] , t + m_c[stage-1]*dt );
+ if( stage == m_s )
+ rt_algebra::foreach( x , x , m_b , m_F , dt );
+ else
+ rt_algebra::foreach( m_x_tmp , x , m_a[stage-1] , m_F , dt );
+ }
+ }
+
+
+private:
+ size_t m_s;
+ vector< vector< double > > m_a;
+ vector< double > m_b;
+ vector< double > m_c;
+
+ state_type m_x_tmp;
+ state_type *m_F;
+};
+
+#endif /* RT_EXPLICIT_RK_HPP_ */