|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r71910 - in sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta: . performance
From: mario.mulansky_at_[hidden]
Date: 2011-05-13 03:31:02
Author: mariomulansky
Date: 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
New Revision: 71910
URL: http://svn.boost.org/trac/boost/changeset/71910
Log:
reorganized rk54ck tests
Added:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk54ck_lorenz.cpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk54ck_lorenz.cpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/lorenz.hpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk54ck_lorenz.cpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk54ck_lorenz.cpp (contents, props changed)
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk54ck_lorenz_def_alg.cpp (contents, props changed)
Removed:
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/gsl_rk4.cpp
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/odeint_rk4.cpp
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/rt_generic_rk4.cpp
Text files modified:
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_explicit_error_rk.hpp | 2
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/Jamfile | 53 +++++++++++++++++++++++++++------------
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk4_lorenz.cpp | 17 +++---------
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk54ck.cpp | 16 ++---------
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk4_lorenz.cpp | 15 +---------
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk54ck.cpp | 16 +----------
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk4_lorenz.cpp | 14 +---------
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk54ck.cpp | 17 +++---------
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_lorenz.cpp | 26 +++----------------
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_lorenz_def_alg.cpp | 26 ++-----------------
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk54ck.cpp | 15 ++---------
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk54ck_def_alg.cpp | 15 ++---------
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/perf_tests.py | 10 +++---
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/results.dat | 2
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/rt_generic_rk4_lorenz.cpp | 15 ++---------
sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/rt_explicit_rk.hpp | 2
16 files changed, 78 insertions(+), 183 deletions(-)
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_explicit_error_rk.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_explicit_error_rk.hpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_explicit_error_rk.hpp 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
@@ -43,7 +43,7 @@
{ }
template< class System >
- void inline do_step( System &system , state_type &x , const double t , const double dt , state_type &x_err )
+ void inline do_step( System system , state_type &x , const double t , const double dt , state_type &x_err )
{
base::do_step( system , x , t , dt );
// compute error estimate
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-13 03:30:59 EDT (Fri, 13 May 2011)
@@ -17,20 +17,20 @@
lib libgsl : : <name>gsl ;
lib libgslcblas : : <name>gslcblas ;
-exe generic_rk4
- : generic_rk4.cpp
- : <toolset>intel:<cxxflags>-inline-forceinline
- ;
+#exe generic_rk4
+# : generic_rk4.cpp
+# : <toolset>intel:<cxxflags>-inline-forceinline
+# ;
exe generic_rk4_lorenz
: generic_rk4_lorenz.cpp
: <toolset>intel:<cxxflags>-inline-forceinline
;
-exe odeint_rk4
- : odeint_rk4.cpp
+#exe odeint_rk4
+# : odeint_rk4.cpp
# : <cxxflags>-Winline
- ;
+# ;
exe odeint_rk4_lorenz
: odeint_rk4_lorenz.cpp
@@ -40,27 +40,27 @@
: odeint_rk4_lorenz_def_alg.cpp
;
-exe nr_rk4
- : nr_rk4.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
+# : 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
+# : gsl_rk4.cpp libgsl libgslcblas
+# ;
exe gsl_rk4_lorenz
: gsl_rk4_lorenz.cpp libgsl libgslcblas
@@ -71,18 +71,39 @@
: <toolset>intel:<cxxflags>-inline-forceinline
;
+exe generic_rk54ck_lorenz
+ : generic_rk54ck_lorenz.cpp
+ : <toolset>intel:<cxxflags>-inline-forceinline
+ ;
+
exe odeint_rk54ck
: odeint_rk54ck.cpp
;
+
+exe odeint_rk54ck_lorenz
+ : odeint_rk54ck_lorenz.cpp
+ ;
exe odeint_rk54ck_def_alg
: odeint_rk54ck_def_alg.cpp
;
+
+exe odeint_rk54ck_lorenz_def_alg
+ : odeint_rk54ck_lorenz_def_alg.cpp
+ ;
exe nr_rk54ck
: nr_rk54ck.cpp
;
+exe nr_rk54ck_lorenz
+ : nr_rk54ck_lorenz.cpp
+ ;
+
exe gsl_rk54ck
: gsl_rk54ck.cpp libgsl libgslcblas
+ ;
+
+exe gsl_rk54ck_lorenz
+ : gsl_rk54ck_lorenz.cpp libgsl libgslcblas
;
\ No newline at end of file
Deleted: 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 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
+++ (empty file)
@@ -1,98 +0,0 @@
-/*
- * 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;
-
-struct lorenz
-{
- 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;
-
-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 );
-
- 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 }};
- //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;
- }
- cout << acc << endl;
- return 0;
-}
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk4_lorenz.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk4_lorenz.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk4_lorenz.cpp 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
@@ -4,22 +4,11 @@
#include "rk_performance_test_case.hpp"
+#include "lorenz.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;
@@ -66,6 +55,8 @@
int main()
{
+ srand( 12312354 );
+
fusion_wrapper stepper;
run( stepper );
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk54ck.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk54ck.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk54ck.cpp 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
@@ -16,6 +16,8 @@
#include "../fusion_explicit_error_rk.hpp"
+#include "lorenz.hpp"
+
#define tab "\t"
using namespace std;
@@ -39,18 +41,6 @@
typedef explicit_error_rk< state_type , 6 > rk54ck_fusion_type;
//typedef explicit_rk< state_type , 6 > rk54ck_fusion_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 )
@@ -92,7 +82,7 @@
timer.restart();
for( size_t i=0 ; i<num_of_steps ; ++i, t+=dt )
{
- rk54ck.do_step( lorenz , x , t , dt , x_err );
+ rk54ck.do_step( lorenz() , x , t , dt , x_err );
if( i % 1000 == 0 ) // simulated stepsize control
dt += (dt*1E-3*rand())/RAND_MAX - dt*5E-4;
}
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk54ck_lorenz.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk54ck_lorenz.cpp 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
@@ -0,0 +1,76 @@
+/*
+ * generic_rk54ck_lorenz.cpp
+ *
+ * Created on: May 12, 2011
+ * Author: mario
+ */
+
+#include <boost/array.hpp>
+
+#include "../fusion_explicit_error_rk.hpp"
+
+#include "rk_performance_test_case.hpp"
+
+#include "lorenz.hpp"
+
+typedef boost::array< double , 3 > state_type;
+typedef explicit_error_rk< state_type , 6 > rk54ck_fusion_type;
+
+
+typedef rk54ck_fusion_type::coef_a_type coef_a_type;
+typedef rk54ck_fusion_type::coef_b_type coef_b_type;
+typedef rk54ck_fusion_type::coef_c_type coef_c_type;
+
+const boost::array< double , 1 > a1 = {{ 0.2 }};
+const boost::array< double , 2 > a2 = {{ 3.0/40.0 , 9.0/40 }};
+const boost::array< double , 3 > a3 = {{ 0.3 , -0.9 , 1.2 }};
+const boost::array< double , 4 > a4 = {{ -11.0/54.0 , 2.5 , -70.0/27.0 , 35.0/27.0 }};
+const boost::array< double , 5 > a5 = {{ 1631.0/55296.0 , 175.0/512.0 , 575.0/13824.0 , 44275.0/110592.0 , 253.0/4096.0 }};
+
+const coef_a_type a = fusion::make_vector( a1 , a2 , a3 , a4 , a5 );
+const coef_b_type b = {{ 37.0/378.0 , 0.0 , 250.0/621.0 , 125.0/594.0 , 0.0 , 512.0/1771.0 }};
+const coef_b_type b2 = {{ b[0]-2825.0/27648.0 , b[1]-0.0 , b[2]-18575.0/48384.0 , b[3]-13525.0/55296.0 , b[4]-277.0/14336.0 , b[5]-0.25 }};
+
+const coef_c_type c = {{ 0.0 , 0.2 , 0.3 , 0.6 , 1.0 , 7.0/8 }};
+
+class fusion_wrapper
+{
+public:
+
+ fusion_wrapper() : m_stepper( a , b , b2 , 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_x_err );
+ //m_t += dt;
+ }
+
+ double state( const size_t i ) const
+ { return m_x[i]; }
+
+private:
+ state_type m_x;
+ state_type m_x_err;
+ double m_t;
+ rk54ck_fusion_type m_stepper;
+};
+
+
+
+int main()
+{
+ srand( 12312354 );
+
+ fusion_wrapper stepper;
+
+ run( stepper );
+}
Deleted: 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 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
+++ (empty file)
@@ -1,87 +0,0 @@
-/*
- * 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( 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;
-}
-
-
-const size_t loops = 20;
-
-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 };
-
- for( size_t n=0 ; n<loops ; ++n )
- {
- 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;
- }
- cout << acc << endl;
- gsl_odeiv_step_free( s );
- return 0;
-
-}
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk4_lorenz.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk4_lorenz.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk4_lorenz.cpp 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
@@ -5,24 +5,13 @@
* Author: mario
*/
-#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
#include "rk_performance_test_case.hpp"
-const size_t dim = 3;
+#include "lorenz.hpp"
-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;
-}
+const size_t dim = 3;
class gsl_wrapper
{
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk54ck.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk54ck.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk54ck.cpp 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
@@ -14,6 +14,8 @@
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
+#include "lorenz.hpp"
+
#define tab "\t"
using namespace std;
@@ -32,20 +34,6 @@
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;
-}
-
-
const size_t loops = 20;
int main()
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk54ck_lorenz.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/gsl_rk54ck_lorenz.cpp 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
@@ -0,0 +1,66 @@
+/*
+ * gsl_rk54ck_lorenz.cpp
+ *
+ * Created on: May 12, 2011
+ * Author: mario
+ */
+
+#include <gsl/gsl_odeiv.h>
+
+#include "rk_performance_test_case.hpp"
+
+#include "lorenz.hpp"
+
+const size_t dim = 3;
+
+class gsl_wrapper
+{
+public:
+
+ gsl_wrapper()
+ {
+ m_s = gsl_odeiv_step_alloc( gsl_odeiv_step_rkck , 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 );
+}
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/lorenz.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/lorenz.hpp 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
@@ -0,0 +1,52 @@
+/*
+ * lorenz.hpp
+ *
+ * Created on: May 12, 2011
+ * Author: mario
+ */
+
+#ifndef LORENZ_HPP_
+#define LORENZ_HPP_
+
+struct lorenz
+{
+ template< class state_type >
+ void inline 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 state_type >
+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];
+}
+
+#include <gsl/gsl_matrix.h>
+
+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;
+}
+
+
+#endif /* LORENZ_HPP_ */
Deleted: 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 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
+++ (empty file)
@@ -1,113 +0,0 @@
-/*
- * 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;
-
-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 )
-{ // 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] );
-}
-
-
-const size_t loops = 20;
-
-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 );
-
- 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_step( lorenz() , x , t , dt );
- acc( timer.elapsed() );
-
- clog.precision( 3 );
- clog.width( 5 );
- clog << acc << " " << x[0] << endl;
- }
- cout << acc << endl;
- return 0;
-}
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk4_lorenz.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk4_lorenz.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk4_lorenz.cpp 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
@@ -9,22 +9,12 @@
#include "rk_performance_test_case.hpp"
+#include "lorenz.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 )
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk54ck.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk54ck.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk54ck.cpp 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
@@ -15,6 +15,8 @@
#include <boost/accumulators/statistics.hpp>
#include <boost/timer.hpp>
+#include "lorenz.hpp"
+
#define tab "\t"
using namespace std;
@@ -33,22 +35,11 @@
typedef boost::timer timer_type;
-
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];
-}
-
template< class System , typename T , size_t dim >
-void rk54ck_step( System &sys , boost::array< T , dim > &x , const double t , const double dt , boost::array< T , dim > &xerr )
+void rk54ck_step( System sys , boost::array< T , dim > &x , const double t , const double dt , boost::array< T , dim > &xerr )
{ // fast rk54ck implementation adapted from the book 'Numerical Recipes'
size_t i;
static const double a2 = 0.2 , a3 = 0.3 , a4 = 0.6 , a5 = 1.0 , a6 = 0.875 ,
@@ -115,7 +106,7 @@
timer.restart();
for( size_t i=0 ; i<num_of_steps ; ++i, t+=dt )
{
- rk54ck_step( lorenz , x , t , dt , x_err );
+ rk54ck_step( lorenz() , x , t , dt , x_err );
if( i % 1000 == 0 ) // simulated stepsize control
dt += (dt*1E-3*rand())/RAND_MAX - dt*5E-4;
}
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk54ck_lorenz.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/nr_rk54ck_lorenz.cpp 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
@@ -0,0 +1,96 @@
+/*
+ * nr_rk54ck_lorenz.cpp
+ *
+ * Created on: May 12, 2011
+ * Author: mario
+ */
+
+#include <boost/array.hpp>
+
+#include "rk_performance_test_case.hpp"
+
+#include "lorenz.hpp"
+
+const size_t dim = 3;
+
+typedef boost::array< double , dim > state_type;
+
+
+template< class System , typename T , size_t dim >
+void rk54ck_step( System sys , boost::array< T , dim > &x , const double t , const double dt , boost::array< T , dim > &xerr )
+{ // fast rk54ck implementation adapted from the book 'Numerical Recipes'
+ size_t i;
+ static const double a2 = 0.2 , a3 = 0.3 , a4 = 0.6 , a5 = 1.0 , a6 = 0.875 ,
+ b21 = 0.2 , b31 = 3.0 / 40.0 , b32 = 9.0 / 40.0 , b41 = 0.3 , b42 =
+ -0.9 , b43 = 1.2 , b51 = -11.0 / 54.0 , b52 = 2.5 , b53 =
+ -70.0 / 27.0 , b54 = 35.0 / 27.0 , b61 = 1631.0 / 55296.0 ,
+ b62 = 175.0 / 512.0 , b63 = 575.0 / 13824.0 , b64 = 44275.0
+ / 110592.0 , b65 = 253.0 / 4096.0 , c1 = 37.0 / 378.0 , c3 =
+ 250.0 / 621.0 , c4 = 125.0 / 594.0 , c6 = 512.0 / 1771.0 ,
+ dc1 = c1 - 2825.0 / 27648.0 , dc3 = c3 - 18575.0 / 48384.0 , dc4 =
+ c4 - 13525.0 / 55296.0 , dc5 = -277.00 / 14336.0 , dc6 = c6
+ - 0.25;
+ const size_t n = dim;
+ boost::array< T , dim > dydx , ak2 , ak3 , ak4 , ak5 , ak6 , ytemp;
+
+ sys( x , dydx , t );
+ for (i=0;i<n;i++)
+ ytemp[i] = x[i] + b21 * dt * dydx[i];
+
+ sys( ytemp , ak2 , t+a2*dt );
+ for (i=0;i<n;i++)
+ ytemp[i] = x[i] + dt*(b31*dydx[i]+b32*ak2[i]);
+
+ sys( ytemp , ak3 , t+a3*dt );
+ for (i=0;i<n;i++)
+ ytemp[i] = x[i] + dt*(b41*dydx[i]+b42*ak2[i]+b43*ak3[i]);
+
+ sys( ytemp , ak4 , t+a4*dt );
+ for (i=0;i<n;i++)
+ ytemp[i] = x[i] + dt*(b51*dydx[i]+b52*ak2[i]+b53*ak3[i]+b54*ak4[i]);
+
+ sys( ytemp, ak5 , t+a5*dt );
+ for (i=0;i<n;i++)
+ ytemp[i] = x[i] + dt*(b61*dydx[i]+b62*ak2[i]+b63*ak3[i]+b64*ak4[i]+b65*ak5[i]);
+
+ sys( ytemp , ak6 , t+a6*dt );
+ for (i=0;i<n;i++)
+ x[i] += dt*(c1*dydx[i]+c3*ak3[i]+c4*ak4[i]+c6*ak6[i]);
+ for (i=0;i<n;i++)
+ xerr[i] = dt*(dc1*dydx[i]+dc3*ak3[i]+dc4*ak4[i]+dc5*ak5[i]+dc6*ak6[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 )
+ {
+ rk54ck_step( lorenz() , m_x , m_t , dt , m_x_err );
+ }
+
+ double state( const size_t i ) const
+ { return m_x[i]; }
+
+private:
+ state_type m_x;
+ state_type m_x_err;
+ double m_t;
+};
+
+
+
+int main()
+{
+ nr_wrapper stepper;
+
+ run( stepper );
+}
Deleted: 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 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
+++ (empty file)
@@ -1,87 +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 , double , state_type , double ,
- boost::numeric::odeint::array_algebra > rk4_odeint_type;
-
-
-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 )
-{
- 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;
-}
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_lorenz.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_lorenz.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_lorenz.cpp 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
@@ -12,32 +12,12 @@
#include "rk_performance_test_case.hpp"
+#include "lorenz.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
{
@@ -69,6 +49,8 @@
int main()
{
+ srand( 12312354 );
+
odeint_wrapper stepper;
run( stepper );
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_lorenz_def_alg.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_lorenz_def_alg.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk4_lorenz_def_alg.cpp 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
@@ -12,31 +12,11 @@
#include "rk_performance_test_case.hpp"
+#include "lorenz.hpp"
+
typedef boost::array< double , 3 > state_type;
typedef boost::numeric::odeint::explicit_rk4< state_type > 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
{
@@ -51,7 +31,7 @@
inline void do_step( const double dt )
{
- m_stepper.do_step( lorenz_func , m_x , m_t , dt );
+ m_stepper.do_step( lorenz() , m_x , m_t , dt );
//m_t += dt;
}
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk54ck.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk54ck.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk54ck.cpp 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
@@ -16,6 +16,8 @@
#include <boost/accumulators/statistics.hpp>
#include <boost/timer.hpp>
+#include "lorenz.hpp"
+
#define tab "\t"
using namespace std;
@@ -41,17 +43,6 @@
boost::numeric::odeint::array_algebra > rk54_ck_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 )
@@ -77,7 +68,7 @@
timer.restart();
for( size_t i=0 ; i<num_of_steps ; ++i, t+=dt )
{
- rk54_ck_odeint.do_step( lorenz , x , t , dt , x_err );
+ rk54_ck_odeint.do_step( lorenz() , x , t , dt , x_err );
if( i % 1000 == 0 ) // simulated stepsize control
dt += (dt*1E-3*rand())/RAND_MAX - dt*5E-4;
}
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk54ck_def_alg.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk54ck_def_alg.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk54ck_def_alg.cpp 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
@@ -16,6 +16,8 @@
#include <boost/accumulators/statistics.hpp>
#include <boost/timer.hpp>
+#include "lorenz.hpp"
+
#define tab "\t"
using namespace std;
@@ -39,17 +41,6 @@
typedef boost::numeric::odeint::explicit_error_rk54_ck< state_type > rk54_ck_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 )
@@ -75,7 +66,7 @@
timer.restart();
for( size_t i=0 ; i<num_of_steps ; ++i, t+=dt )
{
- rk54_ck_odeint.do_step( lorenz , x , t , dt , x_err );
+ rk54_ck_odeint.do_step( lorenz() , x , t , dt , x_err );
if( i % 1000 == 0 ) // simulated stepsize control
dt += (dt*1E-3*rand())/RAND_MAX - dt*5E-4;
}
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk54ck_lorenz.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk54ck_lorenz.cpp 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
@@ -0,0 +1,57 @@
+/*
+ * odeint_rk54ck_lorenz.cpp
+ *
+ * Created on: May 12, 2011
+ * Author: mario
+ */
+
+#include <boost/array.hpp>
+
+#include <boost/numeric/odeint/stepper/explicit_error_rk54_ck.hpp>
+#include <boost/numeric/odeint/algebra/array_algebra.hpp>
+
+#include "rk_performance_test_case.hpp"
+
+#include "lorenz.hpp"
+
+typedef boost::array< double , 3 > state_type;
+typedef boost::numeric::odeint::explicit_error_rk54_ck< state_type , double , state_type , double ,
+ boost::numeric::odeint::array_algebra > rk54_odeint_type;
+
+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_x_err );
+ //m_t += dt;
+ }
+
+ double state( const size_t i ) const
+ { return m_x[i]; }
+
+private:
+ state_type m_x;
+ state_type m_x_err;
+ double m_t;
+ rk54_odeint_type m_stepper;
+};
+
+
+
+int main()
+{
+ srand( 12312354 );
+
+ odeint_wrapper stepper;
+
+ run( stepper );
+}
Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk54ck_lorenz_def_alg.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/odeint_rk54ck_lorenz_def_alg.cpp 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
@@ -0,0 +1,56 @@
+/*
+ * odeint_rk54ck_lorenz.cpp
+ *
+ * Created on: May 12, 2011
+ * Author: mario
+ */
+
+#include <boost/array.hpp>
+
+#include <boost/numeric/odeint/stepper/explicit_error_rk54_ck.hpp>
+#include <boost/numeric/odeint/algebra/array_algebra.hpp>
+
+#include "rk_performance_test_case.hpp"
+
+#include "lorenz.hpp"
+
+typedef boost::array< double , 3 > state_type;
+typedef boost::numeric::odeint::explicit_error_rk54_ck< state_type > rk54_odeint_type;
+
+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_x_err );
+ //m_t += dt;
+ }
+
+ double state( const size_t i ) const
+ { return m_x[i]; }
+
+private:
+ state_type m_x;
+ state_type m_x_err;
+ double m_t;
+ rk54_odeint_type m_stepper;
+};
+
+
+
+int main()
+{
+ srand( 12312354 );
+
+ odeint_wrapper stepper;
+
+ run( stepper );
+}
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-13 03:30:59 EDT (Fri, 13 May 2011)
@@ -2,16 +2,16 @@
from os import system
from os.path import isfile
-toolset = "intel"
+toolset = "gcc-4.5"
-#bin_path = "bin/gcc-4.3/release/"
-bin_path = "bin/intel-linux/release/"
+bin_path = "bin/gcc-4.5/release/"
+#bin_path = "bin/intel-linux/release/"
#bin_path = "bin\\msvc-9.0express\\release\\threading-multi\\"
#extension = ".exe"
extension = ""
-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" ]
+bins = [ "odeint_rk4_lorenz" , "odeint_rk4_lorenz_def_alg" , "generic_rk4_lorenz" , "nr_rk4_lorenz" , "gsl_rk4_lorenz" , "rt_generic_rk4_lorenz" ,
+ "odeint_rk54ck_lorenz" , "odeint_rk54ck_lorenz_def_alg" , "generic_rk54ck_lorenz" , "nr_rk54ck_lorenz" , "gsl_rk54ck_lorenz" ]
results = []
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-13 03:30:59 EDT (Fri, 13 May 2011)
@@ -31,7 +31,7 @@
gcc 4.6.0 | 0.95 | 0.95 | 1.02 | 1.10 | 1.94 | Core i7 870 @ 2.93 GHz
icc 12.0.2 | 1.53 | 1.82 | 1.63 | 1.82 | 1.99 | Core i7 870 @ 2.93 GHz
icc 11.1 | 1.91 | 1.93 | 1.84 | 1.15 | 2.16 | Xeon X5650 @ 2.67 GHz
-gcc 4.4.1 | | | | | | Core2Quad Q6600 @ 2.40GHz
+gcc 4.4.1 | | | | | | Core2Quad Q6600 @ 2.40GHz
msvc 9.0 | 12.8 | | 15.4 | | ---- | Via Nano @ 1.60 GHz
msvc 9.0 | 4.21015 | 4.5484 | 4.79375 | 2.8133 | ---- | Core 2 Quad Q6600 @ 2.40GHz (32bit)
icc 11.1 | 1.80 | 2.44 | 2.39 | | 2.69 | PhenomII X4 945 @ 3 GHz
Deleted: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/rt_generic_rk4.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/rt_generic_rk4.cpp 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
+++ (empty file)
@@ -1,146 +0,0 @@
-/*
- * rt_generic_rk4.cpp
- *
- * Created on: Apr 29, 2011
- * Author: mario
- */
-
-#include <iostream>
-#include <fstream>
-#include <vector>
-
-#include <boost/array.hpp>
-
-#include <boost/accumulators/accumulators.hpp>
-#include <boost/accumulators/statistics.hpp>
-#include <boost/timer.hpp>
-
-#include "rt_algebra.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;
-
-
-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( 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() { 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;
-};
-
-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];
-}
-
-typedef rt_explicit_rk< state_type > rk_stepper_type;
-
-const size_t stage_count = 4;
-
-const size_t loops = 20;
-
-int main( int argc , char **argv )
-{
- 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;
-
- rk_stepper_type rk4_rt( stage_count , a , b , c );
-
- 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 }};
- //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_rt.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;
-}
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/rt_generic_rk4_lorenz.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/rt_generic_rk4_lorenz.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/rt_generic_rk4_lorenz.cpp 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
@@ -11,23 +11,14 @@
#include "rk_performance_test_case.hpp"
+#include "lorenz.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
{
@@ -59,7 +50,7 @@
inline void do_step( const double dt )
{
- m_stepper.do_step( lorenz , m_x , m_t , dt );
+ m_stepper.do_step( lorenz() , m_x , m_t , dt );
//m_t += dt;
}
Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/rt_explicit_rk.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/rt_explicit_rk.hpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/rt_explicit_rk.hpp 2011-05-13 03:30:59 EDT (Fri, 13 May 2011)
@@ -46,7 +46,7 @@
}
template< class System >
- void do_step( System &sys , state_type &x , const double t , const double dt )
+ 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 );
Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk