[Boost-commit] svn:boost r71501 - in sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta: . bin bin/gcc-4.5 bin/gcc-4.5/debug bin/gcc-4.5/release bin/intel-linux bin/intel-linux/release

Subject: [Boost-commit] svn:boost r71501 - in sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta: . bin bin/gcc-4.5 bin/gcc-4.5/debug bin/gcc-4.5/release bin/intel-linux bin/intel-linux/release
From: mario.mulansky_at_[hidden]
Date: 2011-04-26 08:36:37

Author: mariomulansky
Date: 2011-04-26 08:36:34 EDT (Tue, 26 Apr 2011)
New Revision: 71501
URL: http://svn.boost.org/trac/boost/changeset/71501

Log:
fusion runge kutta stepper reorganized in new folder
Added:
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/Jamfile   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/bin/
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/bin/gcc-4.5/
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/bin/gcc-4.5/debug/
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/bin/gcc-4.5/debug/test_explicit_rk   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/bin/gcc-4.5/release/
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/bin/gcc-4.5/release/test_explicit_rk   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/bin/intel-linux/
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/bin/intel-linux/release/
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/bin/intel-linux/release/test_explicit_rk   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_algebra.hpp   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_explicit_rk.hpp   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/test_explicit_rk.cpp   (contents, props changed)

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/Jamfile
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/Jamfile	2011-04-26 08:36:34 EDT (Tue, 26 Apr 2011)
@@ -0,0 +1,20 @@
+# (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 HOME : [ os.environ HOME ] ;
+#path-constant CHRONO_ROOT : [ os.environ CHRONO_ROOT ] ;
+
+project
+    : requirements
+      <define>BOOST_ALL_NO_LIB=1
+      <include>../../../../..
+    ;
+
+exe test_explicit_rk
+    : test_explicit_rk.cpp
+    ;
\ No newline at end of file

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/bin/gcc-4.5/debug/test_explicit_rk
==============================================================================
Binary file. No diff available.

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/bin/gcc-4.5/release/test_explicit_rk
==============================================================================
Binary file. No diff available.

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/bin/intel-linux/release/test_explicit_rk
==============================================================================
Binary file. No diff available.

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_algebra.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_algebra.hpp	2011-04-26 08:36:34 EDT (Tue, 26 Apr 2011)
@@ -0,0 +1,105 @@
+/*
+ * fusion_algebra.hpp
+ *
+ *  Created on: Apr 26, 2011
+ *      Author: mario
+ */
+
+#ifndef FUSION_ALGEBRA_HPP_
+#define FUSION_ALGEBRA_HPP_
+
+#include <boost/array.hpp>
+
+#include <boost/numeric/odeint/algebra/range_algebra.hpp>
+#include <boost/numeric/odeint/algebra/default_operations.hpp>
+
+template< size_t n >
+struct fusion_algebra
+{
+
+    template< class state_type >
+    inline static void foreach( state_type &x_tmp , const state_type &x , const boost::array< double , n > &a ,
+            const state_type *k_vector , const double dt )
+    {
+        for( size_t i=0 ; i<x.size() ; ++i )
+        {
+            x_tmp[i] = x[i];
+            for( size_t j = 0 ; j<n ; ++j )
+                x_tmp[i] += a[j]*dt*k_vector[j][i];
+        }
+    }
+
+};
+
+/*
+
+template<>
+struct fusion_algebra< 1 >
+{
+
+    template< class state_type >
+    inline static void foreach( state_type &x_tmp , const state_type &x , const boost::array< double , 1 > &a ,
+            const state_type *k_vector , const double dt )
+    {
+        for( size_t i=0 ; i<x.size() ; ++i )
+        {
+            x_tmp[i] = x[i] + a[0]*dt*k_vector[0][i];
+        }
+    }
+
+};
+
+
+template<>
+struct fusion_algebra< 2 >
+{
+
+    template< class state_type >
+    inline static void foreach( state_type &x_tmp , const state_type &x , const boost::array< double , 2 > &a ,
+            const state_type *k_vector , const double dt )
+    {
+        for( size_t i=0 ; i<x.size() ; ++i )
+        {
+            x_tmp[i] = x[i] + a[0]*dt*k_vector[0][i] + a[1]*dt*k_vector[1][i];
+        }
+    }
+
+};
+
+
+template<>
+struct fusion_algebra< 3 >
+{
+
+    template< class state_type >
+    inline static void foreach( state_type &x_tmp , const state_type &x , const boost::array< double , 3 > &a ,
+            const state_type *k_vector , const double dt )
+    {
+        for( size_t i=0 ; i<x.size() ; ++i )
+        {
+            x_tmp[i] = x[i] + a[0]*dt*k_vector[0][i] + a[1]*dt*k_vector[1][i] + a[2]*dt*k_vector[2][i];
+        }
+    }
+
+};
+
+
+template<>
+struct fusion_algebra< 4 >
+{
+
+    template< class state_type >
+    inline static void foreach( state_type &x_tmp , const state_type &x , const boost::array< double , 4 > &a ,
+            const state_type *k_vector , const double dt )
+    {
+        for( size_t i=0 ; i<x.size() ; ++i )
+        {
+            x_tmp[i] = x[i] + a[0]*dt*k_vector[0][i] + a[1]*dt*k_vector[1][i] +
+                           a[2]*dt*k_vector[2][i] + a[3]*dt*k_vector[3][i];;
+        }
+    }
+
+};
+*/
+
+#endif /* FUSION_ALGEBRA_HPP_ */

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_explicit_rk.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/fusion_explicit_rk.hpp	2011-04-26 08:36:34 EDT (Tue, 26 Apr 2011)
@@ -0,0 +1,212 @@
+/*
+ * fusion_runge_kutta.hpp
+ *
+ *  Created on: Apr 26, 2011
+ *      Author: mario
+ */
+
+#ifndef FUSION_EXPLICIT_RK_HPP_
+#define FUSION_EXPLICIT_RK_HPP_
+
+#include <boost/mpl/vector.hpp>
+#include <boost/mpl/push_back.hpp>
+#include <boost/mpl/for_each.hpp>
+#include <boost/mpl/range_c.hpp>
+#include <boost/mpl/copy.hpp>
+#include <boost/mpl/size_t.hpp>
+
+#include <boost/fusion/container.hpp>
+#include <boost/fusion/algorithm.hpp>
+#include <boost/fusion/sequence.hpp>
+#include <boost/fusion/adapted.hpp>
+
+
+#include <algorithm>
+#include <iostream>
+#include <string>
+
+#include <boost/array.hpp>
+#include <typeinfo>
+
+#include "fusion_algebra.hpp"
+
+namespace mpl = boost::mpl;
+namespace fusion = boost::fusion;
+
+using namespace std;
+
+struct intermediate_stage {};
+struct last_stage {};
+
+
+
+template< class T , class Constant >
+struct array_wrapper
+{
+    typedef typename boost::array< T , Constant::value > type;
+};
+
+template< class T , class Constant , class StageCategory >
+struct stage_fusion_wrapper
+{
+    typedef typename fusion::vector< size_t , T , boost::array< T , Constant::value > , StageCategory > type;
+};
+
+
+
+template< class StateType , size_t stage_count >
+class explicit_rk
+{
+
+public:
+
+    typedef StateType state_type;
+
+    typedef mpl::range_c< size_t , 1 , stage_count > stage_indices;
+
+    typedef typename fusion::result_of::as_vector
+    <
+        typename mpl::copy
+        <
+            stage_indices ,
+            mpl::inserter
+            <
+                mpl::vector0< > ,
+                mpl::push_back< mpl::_1 , array_wrapper< double , mpl::_2 > >
+            >
+        >::type
+    >::type coef_a_type;
+
+    typedef boost::array< double , stage_count > coef_b_type;
+    typedef boost::array< double , stage_count > coef_c_type;
+
+    typedef typename fusion::result_of::as_vector
+    <
+        typename mpl::push_back
+        <
+            typename mpl::copy
+            <
+                stage_indices,
+                mpl::inserter
+                <
+                    mpl::vector0<> ,
+                    mpl::push_back< mpl::_1 , stage_fusion_wrapper< double , mpl::_2 , intermediate_stage > >
+                >
+            >::type ,
+            typename stage_fusion_wrapper< double , mpl::size_t< stage_count > , last_stage >::type
+        >::type
+    >::type stage_vector_base;
+
+
+    struct stage_vector : public stage_vector_base
+    {
+        struct do_insertion
+        {
+            stage_vector_base &m_base;
+            const coef_a_type &m_a;
+            const coef_c_type &m_c;
+
+            do_insertion( stage_vector_base &base , const coef_a_type &a , const coef_c_type &c )
+            : m_base( base ) , m_a( a ) , m_c( c ) { }
+
+            template< class Index >
+            void operator()( Index ) const
+            {
+                fusion::at_c< 0 >( fusion::at< Index >( m_base ) ) = Index::value;
+                fusion::at_c< 1 >( fusion::at< Index >( m_base ) ) = m_c[ Index::value ];
+                fusion::at_c< 2 >( fusion::at< Index >( m_base ) ) = fusion::at< Index >( m_a );
+            }
+        };
+
+        stage_vector( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c )
+        {
+            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< 1 >( fusion::at_c< stage_count - 1 >( *this ) ) = c[ stage_count - 1 ];
+            fusion::at_c< 2 >( fusion::at_c< stage_count - 1 >( *this ) ) = b;
+        }
+    };
+
+
+
+    template< class System >
+    struct calculate_stage
+    {
+        System &system;
+        state_type &x , &x_tmp;
+        state_type *k_vector;
+        const double t;
+        const double dt;
+
+        calculate_stage( System &_system , state_type &_x , state_type &_x_tmp , state_type *_k_vector ,
+                            const double _t , const double _dt )
+        : system( _system ) , x( _x ) , x_tmp( _x_tmp ) , k_vector( _k_vector ) , t( _t ) , dt( _dt )
+        {}
+
+
+        template< typename T , size_t stage_number >
+        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 );
+
+            if( stage_number == 1 )
+                system( x , k_vector[stage_number-1] , t + c * dt );
+            else
+                system( x_tmp , k_vector[stage_number-1] , t + c * dt );
+
+            fusion_algebra<stage_number>::foreach( x_tmp , x , fusion::at_c< 2 >( stage ) , k_vector , dt);
+        }
+
+
+        template< typename T , size_t stage_number >
+        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 );
+
+            if( stage_number == 1 )
+                system( x , k_vector[stage_number-1] , t + c * dt );
+            else
+                system( x_tmp , k_vector[stage_number-1] , t + c * dt );
+
+            fusion_algebra<stage_number>::foreach( x , x , fusion::at_c< 2 >( stage ) , k_vector , dt);
+        }
+
+
+    };
+
+
+
+
+public:
+
+    explicit_rk( const coef_a_type &a ,
+                  const coef_b_type &b ,
+                  const coef_c_type &c )
+    : m_a( a ) , m_b( b ) , m_c( c ) ,
+      m_stages( a , b , c )
+
+    { }
+
+
+    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 ) );
+    }
+
+
+private:
+
+    const coef_a_type m_a;
+    const coef_b_type m_b;
+    const coef_c_type m_c;
+    const stage_vector m_stages;
+    state_type m_x_tmp;
+    state_type m_k_vector[stage_count];
+};
+
+
+#endif /* FUSION_EXPLICIT_RK_HPP_ */

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/test_explicit_rk.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/test_explicit_rk.cpp	2011-04-26 08:36:34 EDT (Tue, 26 Apr 2011)
@@ -0,0 +1,102 @@
+/*
+ * test_explicit_rk.cpp
+ *
+ *  Created on: Apr 26, 2011
+ *      Author: mario
+ */
+
+#include <iostream>
+#include <fstream>
+
+#include <tr1/array>
+
+#include <boost/accumulators/accumulators.hpp>
+#include <boost/accumulators/statistics.hpp>
+#include <boost/timer.hpp>
+
+#include "fusion_explicit_rk.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 std::tr1::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 size_t dt = 0.01;
+
+    accumulator_type acc;
+    timer_type timer;
+
+    srand48( 12312354 );
+
+    while( true )
+    {
+        state_type x = {{ 10.0 * drand48() , 10.0 * drand48() , 10.0 * drand48() }};
+        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;
+}
+
+/*
+ * Compile with :
+ * g++ -O3 -I$BOOST_ROOT -I$HOME/boost/chrono -I$ODEINT_ROOT butcher_performance.cpp
+ */