Boost logo

Boost-Commit :

Subject: [Boost-commit] svn:boost r72241 - in sandbox/odeint/branches/karsten: boost/numeric/odeint/algebra boost/numeric/odeint/algebra/detail boost/numeric/odeint/stepper libs/numeric/odeint/ideas/fusion_runge_kutta/performance libs/numeric/odeint/test
From: mario.mulansky_at_[hidden]
Date: 2011-05-28 08:24:48


Author: mariomulansky
Date: 2011-05-28 08:24:46 EDT (Sat, 28 May 2011)
New Revision: 72241
URL: http://svn.boost.org/trac/boost/changeset/72241

Log:
improved generic stepper implementation - now uses range algebra. however, code is still a bit ugly...
Text files modified:
   sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/default_operations.hpp | 40 ++++++++++++++++++++++++++++
   sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/detail/for_each.hpp | 10 +++++-
   sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/range_algebra.hpp | 56 ++++++++++++++++++++++++++++++++++++++++
   sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_generic_rk.hpp | 11 ++++---
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/fusion_runge_kutta/performance/generic_rk4_lorenz.cpp | 2
   sandbox/odeint/branches/karsten/libs/numeric/odeint/test/generic_stepper.cpp | 15 +++++++---
   6 files changed, 121 insertions(+), 13 deletions(-)

Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/default_operations.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/default_operations.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/default_operations.hpp 2011-05-28 08:24:46 EDT (Sat, 28 May 2011)
@@ -227,11 +227,51 @@
         };
 
 
+ /*
+ * For generic case
+ */
+
+ template< size_t n , class Fac = double >
+ struct scale_sumn
+ {
+ BOOST_STATIC_ASSERT( false );
+ };
+
+ template< class Fac >
+ struct scale_sumn< 1 , Fac > : public scale_sum2< Fac >
+ {
+ scale_sumn( const boost::array<Fac,1> &a , const Fac &dt ) : scale_sum2( 1.0 , a[0]*dt )
+ { }
+
+ typedef void result_type;
+ };
 
+ template< class Fac >
+ struct scale_sumn< 2 , Fac > : public scale_sum3< Fac >
+ {
+ scale_sumn( const boost::array<Fac,2> &a , const Fac &dt ) : scale_sum3( 1.0 , a[0]*dt , a[1]*dt )
+ { }
 
+ typedef void result_type;
+ };
+
+ template< class Fac >
+ struct scale_sumn< 3 , Fac > : public scale_sum4< Fac >
+ {
+ scale_sumn( const boost::array<Fac,3> &a , const Fac &dt ) : scale_sum4( 1.0 , a[0]*dt , a[1]*dt , a[2]*dt )
+ { }
 
+ typedef void result_type;
+ };
 
+ template< class Fac >
+ struct scale_sumn< 4 , Fac > : public scale_sum5< Fac >
+ {
+ scale_sumn( const boost::array<Fac,4> &a , const Fac &dt ) : scale_sum5( 1.0 , a[0]*dt , a[1]*dt , a[2]*dt , a[3]*dt )
+ { }
 
+ typedef void result_type;
+ };
 
         /*
          * for usage in for_each2

Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/detail/for_each.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/detail/for_each.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/detail/for_each.hpp 2011-05-28 08:24:46 EDT (Sat, 28 May 2011)
@@ -84,8 +84,14 @@
         for( ; first1 != last1 ; )
                 op( *first1++ , *first2++ , *first3++ , *first4++ , *first5++ , *first6++ , *first7++ , *first8++ );
 }
-
-
+//
+//template< size_t n , class Iterator1 , class Iterator2 , class Iterator3 , class Iterator4 , class Operation >
+//inline void for_eachn( Iterator1 first1 , Iterator1 last1 , Iterator2 first2 , Iterator4 first3 , Iterator4 first_array[n-1] , Operation op )
+//{
+// while( first1 != last1 )
+// {
+// op( *first1++ , *first2++ , *first3++ , first_array ,
+//}
 
 } // detail
 } // odeint

Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/range_algebra.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/range_algebra.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/algebra/range_algebra.hpp 2011-05-28 08:24:46 EDT (Sat, 28 May 2011)
@@ -14,6 +14,7 @@
 #define BOOST_BOOST_NUMERIC_ODEINT_STANDARD_ALGEBRA_HPP_INCLUDED
 
 #include <boost/range.hpp>
+#include <boost/mpl/size_t.hpp>
 
 #include <boost/numeric/odeint/algebra/detail/macros.hpp>
 #include <boost/numeric/odeint/algebra/detail/for_each.hpp>
@@ -82,6 +83,61 @@
         }
 
 
+
+ /* for the generic stepper
+ */
+ template< size_t n , class S1 , class S2 , class S3 , class S4 , class Op >
+ inline static void for_eachn( S1 &s1 , S2 &s2 , S3 &s3 , S4 s4_array[n] , Op op )
+ {
+ for_eachn_fw( s1 , s2 , s3 , s4_array , op , boost::mpl::size_t< n-1 >() );
+ }
+
+ template< class S1 , class S2 , class S3 , class S4 , class Op >
+ inline static void for_eachn_fw( S1 &s1 , S2 &s2 , S3 &s3 , S4 s4_array[0] , Op op ,
+ boost::mpl::size_t< 0 > c )
+ {
+ detail::for_each3( boost::begin( s1 ) , boost::end( s1 ) , boost::begin( s2 ) ,
+ boost::begin( s3 ) , op );
+ }
+
+ template< class S1 , class S2 , class S3 , class S4 , class Op >
+ inline static void for_eachn_fw( S1 &s1 , S2 &s2 , S3 &s3 , S4 s4_array[1] , Op op ,
+ boost::mpl::size_t< 1 > c )
+ {
+ detail::for_each4( boost::begin( s1 ) , boost::end( s1 ) , boost::begin( s2 ) ,
+ boost::begin( s3 ) , boost::begin( s4_array[0] ) , op );
+ }
+
+ template< class S1 , class S2 , class S3 , class S4 , class Op >
+ inline static void for_eachn_fw( S1 &s1 , S2 &s2 , S3 &s3 , S4 s4_array[2] , Op op ,
+ boost::mpl::size_t< 2 > c )
+ {
+ detail::for_each5( boost::begin( s1 ) , boost::end( s1 ) , boost::begin( s2 ) ,
+ boost::begin( s3 ) , boost::begin( s4_array[0] ) ,
+ boost::begin( s4_array[1] ) , op );
+ }
+
+ template< class S1 , class S2 , class S3 , class S4 , class Op >
+ inline static void for_eachn_fw( S1 &s1 , S2 &s2 , S3 &s3 , S4 s4_array[3] , Op op ,
+ boost::mpl::size_t< 3 > c )
+ {
+ detail::for_each6( boost::begin( s1 ) , boost::end( s1 ) , boost::begin( s2 ) ,
+ boost::begin( s3 ) , boost::begin( s4_array[0] ) ,
+ boost::begin( s4_array[1] ) , boost::begin( s4_array[2] ) , op );
+ }
+
+ /*template< size_t n , class StateOut , class StateIn , class DerivIn , typename T >
+ inline static void foreach( StateOut &x_tmp ,
+ const boost::array< T , n > &a ,
+ const DerivIn &dxdt , const StateIn k_vector[n] , const T dt )
+ {
+ for( size_t i=0 ; i<x.size() ; ++i )
+ {
+ x_tmp[i] = a[0]*dt*dxdt[i];
+ for( size_t j = 1 ; j<n ; ++j )
+ x_tmp[i] += a[j]*dt*k_vector[j-1][i];
+ }
+ }*/
 };
 
 } // odeint

Modified: sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_generic_rk.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_generic_rk.hpp (original)
+++ sandbox/odeint/branches/karsten/boost/numeric/odeint/stepper/explicit_generic_rk.hpp 2011-05-28 08:24:46 EDT (Sat, 28 May 2011)
@@ -69,7 +69,7 @@
     class Value = double ,
     class Deriv = State ,
     class Time = Value ,
- class Algebra = generic_algebra ,
+ class Algebra = range_algebra ,
         class Operations = default_operations ,
         class AdjustSizePolicy = adjust_size_initially_tag
>
@@ -196,9 +196,11 @@
                         //std::cout << stage_number-2 << ", t': " << t + stage.c * dt << std::endl;
 
                         if( stage_number < StageCount )
- generic_algebra::foreach<stage_number>( x_tmp , x , stage.a , dxdt , F , dt);
+ algebra_type::for_eachn<stage_number>( x_tmp , x , dxdt , F ,
+ typename operations_type::template scale_sumn< stage_number , time_type >( stage.a , dt ) );
                         else
- generic_algebra::foreach<stage_number>( x_out , x , stage.a , dxdt , F , dt);
+ algebra_type::for_eachn<stage_number>( x_out , x , dxdt , F ,
+ typename operations_type::template scale_sumn< stage_number , time_type >( stage.a , dt ) );
         }
 
     };
@@ -233,5 +235,4 @@
 }
 }
 }
-
-#endif /* EXPLICIT_GENERIC_RK_HPP_ */
+#endif /* EXPLICIT_GENERIC_RK_HPP_ */
\ No newline at end of file

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-28 08:24:46 EDT (Sat, 28 May 2011)
@@ -13,7 +13,7 @@
 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 }};
 

Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/test/generic_stepper.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/test/generic_stepper.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/test/generic_stepper.cpp 2011-05-28 08:24:46 EDT (Sat, 28 May 2011)
@@ -17,6 +17,7 @@
 #include <boost/test/unit_test.hpp>
 
 #include <boost/numeric/odeint/stepper/explicit_generic_rk.hpp>
+#include <boost/numeric/odeint/stepper/explicit_rk4.hpp>
 
 #include <boost/array.hpp>
 
@@ -42,28 +43,32 @@
 const stepper_type::coef_b_type b = {{ 1.0/6 , 1.0/3 , 1.0/3 , 1.0/6 }};
 const stepper_type::coef_c_type c = {{ 0.0 , 0.5 , 0.5 , 1.0 }};
 
+typedef explicit_rk4< state_type > rk4_stepper_type;
+
 BOOST_AUTO_TEST_SUITE( generic_stepper_test )
 
 BOOST_AUTO_TEST_CASE( test_generic_stepper )
 {
         stepper_type stepper( a , b , c );
 
+ rk4_stepper_type rk4;
+
         typedef stepper_type::state_type state_type;
         typedef stepper_type::value_type stepper_value_type;
         typedef stepper_type::deriv_type deriv_type;
         typedef stepper_type::time_type time_type;
 
         state_type x = {{ 0.0 , 1.0 }};
+ state_type y = x;
         
         stepper.do_step( sys , x , 0.0 , 0.1 );
 
- /*using std::abs;
- value_type eps = 1E-12;
+ rk4.do_step( sys , y , 0.0 , 0.1 );
 
         // compare with analytic solution of above system
- BOOST_CHECK_SMALL( abs( x[0] - 20.0/81.0 ) , eps );
- BOOST_CHECK_SMALL( abs( x[1] - 10.0/9.0 ) , eps );
- */
+ BOOST_CHECK_EQUAL( x[0] , y[0] );
+ BOOST_CHECK_EQUAL( x[1] , y[1] );
+
 }
 
 BOOST_AUTO_TEST_SUITE_END()
\ No newline at end of file


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