[Boost-commit] svn:boost r71560 - in sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas: fusion_runge_kutta fusion_runge_kutta/performance generic_stepper

Subject: [Boost-commit] svn:boost r71560 - in sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas: fusion_runge_kutta fusion_runge_kutta/performance generic_stepper
From: mario.mulansky_at_[hidden]
Date: 2011-04-28 09:38:16

Author: mariomulansky
Date: 2011-04-28 09:38:14 EDT (Thu, 28 Apr 2011)
New Revision: 71560
URL: http://svn.boost.org/trac/boost/changeset/71560

Log:
reorganized generic fusion stepper
performance tests
Added:
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/Jamfile   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk4.cpp   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk4.cpp   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk4.cpp   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4.cpp   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/results.dat   (contents, props changed)
Text files modified: 
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/Jamfile                    |     2                                         
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_algebra.hpp         |     4 +-                                      
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_explicit_rk_new.hpp |    47 ++++++++++++++++----------------------- 
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/test_explicit_rk.cpp       |     4 +-                                      
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/Jamfile                       |     2                                         
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/fusion_algebra.hpp            |     1                                         
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance.cpp               |     2                                         
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance_fusion.cpp        |     3 -                                       
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance_odeint.cpp        |     7 -----                                   
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/runge_kutta_stepper_fast.hpp  |     8 +++---                                  
   10 files changed, 32 insertions(+), 48 deletions(-)

Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/Jamfile	(original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/Jamfile	2011-04-28 09:38:14 EDT (Thu, 28 Apr 2011)
@@ -7,7 +7,7 @@
 import path ; 
 
 #path-constant HOME : [ os.environ HOME ] ;
-#path-constant CHRONO_ROOT : [ os.environ CHRONO_ROOT ] ;
+path-constant CHRONO_ROOT : [ os.environ CHRONO_ROOT ] ;
 
 project
     : requirements

Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_algebra.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_algebra.hpp	(original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_algebra.hpp	2011-04-28 09:38:14 EDT (Thu, 28 Apr 2011)
@@ -38,7 +38,7 @@
 
 /* !!!!!!!   Actually, this is factor 3 slower with intel compiler, so we don'y use it !!!!!
  * Update: Current implementation increases performance on msvc 9.0 by about 30%, so it is in use again....
- */
+ *
 
 template<>
 struct fusion_algebra< 1 >
@@ -109,6 +109,6 @@
     }
 
 };
-
+*/
 
 #endif /* FUSION_ALGEBRA_HPP_ */

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-04-28 09:38:14 EDT (Thu, 28 Apr 2011)
@@ -55,16 +55,16 @@
 template< class T , size_t i , class StageCategory >
 struct stage
 {
-    T m_c;
-    boost::array< T , i > m_a;
+    T c;
+    boost::array< T , i > a;
     typedef StageCategory category;
 };
 
 template< class T , size_t i>
 struct stage< T , i , last_stage >
 {
-    T m_c;
-    boost::array< T , i > m_b;
+    T c;
+    boost::array< T , i > b;
     typedef last_stage category;
 };
 
@@ -136,8 +136,8 @@
             void operator()( Index ) const
             {
                 //fusion::at< Index >( m_base )::  = Index::value;
-                fusion::at< Index >( m_base ).m_c  = m_c[ Index::value ];
-                fusion::at< Index >( m_base ).m_a = fusion::at< Index >( m_a );
+                fusion::at< Index >( m_base ).c  = m_c[ Index::value ];
+                fusion::at< Index >( m_base ).a = fusion::at< Index >( m_a );
             }
         };
 
@@ -146,8 +146,8 @@
             typedef mpl::range_c< size_t , 0 , stage_count - 1 > indices;
             mpl::for_each< indices >( do_insertion( *this , a , c ) );
             //fusion::at_c< 0 >( fusion::at_c< stage_count - 1 >( *this ) ) = stage_count - 1 ;
-            fusion::at_c< stage_count - 1 >( *this ).m_c = c[ stage_count - 1 ];
-            fusion::at_c< stage_count - 1 >( *this ).m_b = b;
+            fusion::at_c< stage_count - 1 >( *this ).c = c[ stage_count - 1 ];
+            fusion::at_c< stage_count - 1 >( *this ).b = b;
         }
     };
 
@@ -158,13 +158,13 @@
     {
         System &system;
         state_type &x , &x_tmp;
-        state_type *k_vector;
+        state_type *F;
         const double t;
         const double dt;
 
-        calculate_stage( System &_system , state_type &_x , state_type &_x_tmp , state_type *_k_vector ,
+        calculate_stage( System &_system , state_type &_x , state_type &_x_tmp , state_type *_F ,
                             const double _t , const double _dt )
-        : system( _system ) , x( _x ) , x_tmp( _x_tmp ) , k_vector( _k_vector ) , t( _t ) , dt( _dt )
+        : system( _system ) , x( _x ) , x_tmp( _x_tmp ) , F( _F ) , t( _t ) , dt( _dt )
         {}
 
 
@@ -172,14 +172,12 @@
         void operator()( stage< T , stage_number , intermediate_stage > const &stage ) const
         //typename stage_fusion_wrapper< T , mpl::size_t< stage_number > , intermediate_stage >::type const &stage ) const
         {
-            double c = stage.m_c;
-
             if( stage_number == 1 )
-                system( x , k_vector[stage_number-1] , t + c * dt );
+                system( x , F[stage_number-1] , t + stage.c * dt );
             else
-                system( x_tmp , k_vector[stage_number-1] , t + c * dt );
+                system( x_tmp , F[stage_number-1] , t + stage.c * dt );
 
-            fusion_algebra<stage_number>::foreach( x_tmp , x , stage.m_a , k_vector , dt);
+            fusion_algebra<stage_number>::foreach( x_tmp , x , stage.a , F , dt);
         }
 
 
@@ -187,22 +185,17 @@
         void operator()( stage< T , stage_number , last_stage > const &stage ) const
         //void operator()( typename stage_fusion_wrapper< T , mpl::size_t< stage_number > , last_stage >::type const &stage ) const
         {
-            double c = stage.m_c;
-
             if( stage_number == 1 )
-                system( x , k_vector[stage_number-1] , t + c * dt );
+                system( x , F[stage_number-1] , t + stage.c * dt );
             else
-                system( x_tmp , k_vector[stage_number-1] , t + c * dt );
+                system( x_tmp , F[stage_number-1] , t + stage.c * dt );
 
-            fusion_algebra<stage_number>::foreach( x , x , stage.m_b , k_vector , dt);
+            fusion_algebra<stage_number>::foreach( x , x , stage.b , F , dt);
         }
 
 
     };
 
-
-
-
 public:
 
     explicit_rk( const coef_a_type &a ,
@@ -217,10 +210,9 @@
     template< class System >
     void do_step( System &system , state_type &x , double t , const double dt )
     {
-        fusion::for_each( m_stages , calculate_stage< System >( system , x , m_x_tmp , m_k_vector , t , dt ) );
+        fusion::for_each( m_stages , calculate_stage< System >( system , x , m_x_tmp , m_F , t , dt ) );
     }
 
-
 private:
 
     const coef_a_type m_a;
@@ -228,8 +220,7 @@
     const coef_c_type m_c;
     const stage_vector m_stages;
     state_type m_x_tmp;
-    state_type m_k_vector[stage_count];
+    state_type m_F[stage_count];
 };
 
-
 #endif /* FUSION_EXPLICIT_RK_HPP_ */

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/Jamfile
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/Jamfile	2011-04-28 09:38:14 EDT (Thu, 28 Apr 2011)
@@ -0,0 +1,35 @@
+# (C) Copyright 2010 : Karsten Ahnert, 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)
+
+import os ;
+import modules ;
+import path ; 
+
+path-constant CHRONO_ROOT : [ os.environ CHRONO_ROOT ] ;
+
+project
+    : requirements
+      <define>BOOST_ALL_NO_LIB=1
+      <include>../../../../../..
+      <include>$(CHRONO_ROOT)
+    ;
+
+exe generic_rk4 
+    : generic_rk4.cpp
+    ;
+   
+exe odeint_rk4
+    :   odeint_rk4.cpp
+    ;
+    
+exe nr_rk4
+    : nr_rk4.cpp
+    ;
+    
+lib libgsl : : <name>gsl ;
+lib libgslcblas : : <name>gslcblas ;   
+    
+exe gsl_rk4
+    : gsl_rk4.cpp libgsl libgslcblas
+    ;
\ No newline at end of file

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk4.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk4.cpp	2011-04-28 09:38:14 EDT (Thu, 28 Apr 2011)
@@ -0,0 +1,97 @@
+/*
+ * generic_rk4.cpp
+ *
+ *  Created on: Apr 28, 2011
+ *      Author: mario
+ */
+
+#include <iostream>
+#include <fstream>
+
+#include <boost/array.hpp>
+
+#include <boost/accumulators/accumulators.hpp>
+#include <boost/accumulators/statistics.hpp>
+#include <boost/timer.hpp>
+
+//#include "fusion_explicit_rk.hpp"
+#include "../fusion_explicit_rk_new.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 explicit_rk< state_type , 4 > rk4_fusion_type;
+
+
+void lorenz( const state_type &x , state_type &dxdt , 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];
+}
+
+
+
+
+int main( int argc , char **argv )
+{
+    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 }};
+
+    rk4_fusion_type rk4_fusion( a , b , c );
+
+    const size_t num_of_steps = 20000000;
+    const double dt = 1E-10;
+
+    accumulator_type acc;
+    timer_type timer;
+
+    srand( 12312354 );
+
+    while( true )
+    {
+        state_type x = {{ 10.0 * rand()/RAND_MAX , 10.0 * rand()/RAND_MAX , 10.0 * rand()/RAND_MAX }};
+        //state_type x = {{ 10.0 , 1.0 , 5.0 }};
+        double t = 0.0;
+
+        timer.restart();
+        for( size_t i=0 ; i<num_of_steps ; ++i, t+=dt )
+            rk4_fusion.do_step( lorenz , x , t , dt );
+        acc( timer.elapsed() );
+
+        clog.precision( 3 );
+        clog.width( 5 );
+        clog << acc << " " << x[0] << endl;
+    }
+
+    return 0;
+}

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk4.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk4.cpp	2011-04-28 09:38:14 EDT (Thu, 28 Apr 2011)
@@ -0,0 +1,84 @@
+/*
+ * gsl_rk4.cpp
+ *
+ *  Created on: Apr 28, 2011
+ *      Author: mario
+ */
+
+#include <iostream>
+
+#include <boost/accumulators/accumulators.hpp>
+#include <boost/accumulators/statistics.hpp>
+#include <boost/timer.hpp>
+
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_odeiv.h>
+
+#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;
+
+
+int lorenz_gsl( 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;
+}
+
+
+int main()
+{
+    const size_t num_of_steps = 20000000 / 3 ; // gsl rk4 routine makes error control by
+                                               // additional doing two steps with half step size
+    const double dt = 1E-10 * 3 ;             // so it actually does 3 * num_of_steps steps
+
+    accumulator_type acc;
+    timer_type timer;
+
+    srand( 12312354 );
+
+    gsl_odeiv_step *s = gsl_odeiv_step_alloc( gsl_odeiv_step_rk4 , 3);
+    gsl_odeiv_system sys = { lorenz_gsl , 0 , 3 , 0 };
+
+    while( true )
+    {
+        double x[3] = { 10.0 * rand()/RAND_MAX ,
+                         10.0 * rand()/RAND_MAX ,
+                         10.0 * rand()/RAND_MAX };
+        double x_err[3];
+
+        double t = 0.0;
+
+        timer.restart();
+        for( size_t i=0 ; i<num_of_steps ; ++i, t+=dt )
+            gsl_odeiv_step_apply ( s , t , dt , x , x_err , 0 , 0 , &sys );
+        acc( timer.elapsed() );
+
+        clog.precision( 3 );
+        clog.width( 5 );
+        clog << acc << " " << x[0] << endl;
+    }
+
+    return 0;
+
+}

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk4.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk4.cpp	2011-04-28 09:38:14 EDT (Thu, 28 Apr 2011)
@@ -0,0 +1,107 @@
+/*
+ * nr_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/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;
+
+void lorenz( const state_type &x , state_type &dxdt , 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];
+}
+
+
+template< class System , typename T , size_t dim >
+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;
+    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] );
+}
+
+
+int main( int argc , char **argv )
+{
+    const size_t num_of_steps = 20000000;
+    const double dt = 1E-10;
+
+    accumulator_type acc;
+    timer_type timer;
+
+    srand( 12312354 );
+
+    while( true )
+    {
+        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_step( lorenz , x , t , dt );
+        acc( timer.elapsed() );
+
+        clog.precision( 3 );
+        clog.width( 5 );
+        clog << acc << " " << x[0] << endl;
+    }
+
+    return 0;
+}

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4.cpp	2011-04-28 09:38:14 EDT (Thu, 28 Apr 2011)
@@ -0,0 +1,87 @@
+/*
+ * 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;
+typedef boost::numeric::odeint::explicit_rk4< state_type , double , state_type , double ,
+                                              boost::numeric::odeint::array_algebra > rk4_odeint_type;
+
+
+void lorenz( const state_type &x , state_type &dxdt , 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];
+}
+
+
+
+
+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 );
+
+    while( true )
+    {
+        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;
+    }
+
+    return 0;
+}

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/results.dat
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/results.dat	2011-04-28 09:38:14 EDT (Thu, 28 Apr 2011)
@@ -0,0 +1,8 @@
+Results for Runge-Kutta 4 with different implementations and compilers
+
+            |   odeint  |   generic |   nr      |   gsl     |
+-------------------------------------------------------------
+gcc 4.5.0   |   0.79    |   0.80    |   0.82    |   1.11    |   Corei7 870 @ 2.93 GHz
+gcc 4.3.2   |   1.00    |   1.07    |   0.97    |   2.03    |   Core2Quad Q9550 @ 2.83 GHz
+icc 12.0.2  |   0.95    |   0.87    |   1.09    |   1.12    |   Corei7 870 @ 2.93 GHz
+icc 11.1    |   1.11    |   0.98    |   1.19    |   1.28    |   Xeon X5650   @ 2.67 GHz
\ No newline at end of file

Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/test_explicit_rk.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/test_explicit_rk.cpp	(original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/test_explicit_rk.cpp	2011-04-28 09:38:14 EDT (Thu, 28 Apr 2011)
@@ -14,8 +14,8 @@
 #include <boost/accumulators/statistics.hpp>
 #include <boost/timer.hpp>
 
-#include "fusion_explicit_rk.hpp"
-//#include "fusion_explicit_rk_new.hpp"
+//#include "fusion_explicit_rk.hpp"
+#include "fusion_explicit_rk_new.hpp"
 
 #define tab "\t"
 

Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/Jamfile	(original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/Jamfile	2011-04-28 09:38:14 EDT (Thu, 28 Apr 2011)
@@ -6,7 +6,7 @@
 import modules ;
 import path ; 
 
-#path-constant CHRONO_ROOT : [ os.environ CHRONO_ROOT ] ;
+path-constant CHRONO_ROOT : [ os.environ CHRONO_ROOT ] ;
 
 project
     : requirements

Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/fusion_algebra.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/fusion_algebra.hpp	(original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/fusion_algebra.hpp	2011-04-28 09:38:14 EDT (Thu, 28 Apr 2011)
@@ -29,6 +29,7 @@
 
 };
 
+
 template<>
 struct fusion_algebra< 1 >
 {

Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance.cpp	(original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance.cpp	2011-04-28 09:38:14 EDT (Thu, 28 Apr 2011)
@@ -23,7 +23,7 @@
 #include <boost/timer.hpp>
 
 #include "predefined_steppers.hpp"
-#include "runge_kutta_stepper.hpp"
+#include "runge_kutta_stepper_fast.hpp"
 
 #define tab "\t"
 

Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance_fusion.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance_fusion.cpp	(original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance_fusion.cpp	2011-04-28 09:38:14 EDT (Thu, 28 Apr 2011)
@@ -21,7 +21,7 @@
 #include <boost/accumulators/statistics.hpp>
 #include <boost/timer.hpp>
 
-#include "runge_kutta_stepper_fast.hpp"
+#include "runge_kutta_stepper.hpp"
 
 #define tab "\t"
 
@@ -58,7 +58,6 @@
 
 
 
-
 int main( int argc , char **argv )
 {
     typedef rk4_fusion_type::coef_a_type coef_a_type;

Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance_odeint.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance_odeint.cpp	(original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/performance_odeint.cpp	2011-04-28 09:38:14 EDT (Thu, 28 Apr 2011)
@@ -5,13 +5,6 @@
  *      Author: mario
  */
 
-/*
- * butcher_test.cpp
- *
- *  Created on: Nov 5, 2010
- *      Author: karsten
- */
-
 #include <iostream>
 #include <fstream>
 

Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/runge_kutta_stepper_fast.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/runge_kutta_stepper_fast.hpp	(original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/runge_kutta_stepper_fast.hpp	2011-04-28 09:38:14 EDT (Thu, 28 Apr 2011)
@@ -135,11 +135,11 @@
         {}
 
 
-        template< typename T , size_t stage_number >
+        template< typename T , const size_t stage_number >
         inline void operator()( fusion::vector< size_t , T , boost::array< T , stage_number > , intermediate_stage > const &stage ) const
         //typename stage_fusion_wrapper< T , mpl::size_t< stage_number > , intermediate_stage >::type const &stage ) const
         {
-            double c = fusion::at_c< 1 >( stage );
+            const double c = fusion::at_c< 1 >( stage );
 
             if( stage_number == 1 )
                 system( x , k_vector[stage_number-1] , t + c * dt );
@@ -150,11 +150,11 @@
         }
 
 
-        template< typename T , size_t stage_number >
+        template< typename T , const size_t stage_number >
         inline void operator()( fusion::vector< size_t , T , boost::array< T , stage_number > , last_stage > const &stage ) const
         //void operator()( typename stage_fusion_wrapper< T , mpl::size_t< stage_number > , last_stage >::type const &stage ) const
         {
-            double c = fusion::at_c< 1 >( stage );
+            const double c = fusion::at_c< 1 >( stage );
 
             if( stage_number == 1 )
                 system( x , k_vector[stage_number-1] , t + c * dt );