[Boost-commit] svn:boost r57716 - in sandbox/odeint: boost/numeric/odeint libs/numeric/odeint/examples

Subject: [Boost-commit] svn:boost r57716 - in sandbox/odeint: boost/numeric/odeint libs/numeric/odeint/examples
From: mario.mulansky_at_[hidden]
Date: 2009-11-16 17:18:45

Author: mariomulansky
Date: 2009-11-16 17:18:44 EST (Mon, 16 Nov 2009)
New Revision: 57716
URL: http://svn.boost.org/trac/boost/changeset/57716

Log:
+ generic RK with c-styled parameter handling
unfortunately not faster...
Text files modified: 
   sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp             |     2                                         
   sandbox/odeint/boost/numeric/odeint/stepper_rk_generic.hpp                |   186 ++++++++++++++++++++++++++++++++++++++++
   sandbox/odeint/libs/numeric/odeint/examples/lorenz_cmp_rk4_rk_generic.cpp |    21 +++                                     
   3 files changed, 203 insertions(+), 6 deletions(-)

Modified: sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp	(original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_rk4_classical.hpp	2009-11-16 17:18:44 EST (Mon, 16 Nov 2009)
@@ -127,8 +127,6 @@
         }
 
 
-        /* RK4 step with error estimation to 5th order according to Cash Karp */
-
     };
 
 } // namespace odeint

Modified: sandbox/odeint/boost/numeric/odeint/stepper_rk_generic.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_rk_generic.hpp	(original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_rk_generic.hpp	2009-11-16 17:18:44 EST (Mon, 16 Nov 2009)
@@ -238,6 +238,192 @@
 
     };
 
+
+
+    /* ############################################################
+       ############################################################
+
+       C-Array Version of a,b,c handling
+
+       ############################################################
+       ############################################################
+    */
+
+
+
+
+    template<
+        class Container ,
+        class Time = double ,
+        class Resizer = resizer< Container >
+        >
+    class stepper_rk_generic_ptr {
+
+        // provide basic typedefs
+    public:
+
+        typedef Container container_type;
+        typedef Resizer resizer_type;
+        typedef Time time_type;
+        typedef const unsigned short order_type;
+        typedef typename container_type::value_type value_type;
+        typedef typename container_type::iterator iterator;
+
+
+        // check the concept of the ContainerType
+    private:
+
+        BOOST_CLASS_REQUIRE( container_type ,
+                             boost::numeric::odeint, Container );
+
+        // private variables
+    private:
+        typedef std::vector< container_type > container_vector;
+        typedef std::vector< iterator > container_iterator_vector;
+
+        container_vector m_xvec;
+        container_iterator_vector m_xiter_vec;
+        container_type m_xtmp;
+        const time_type* m_a;
+        const time_type* m_b;
+        const time_type* m_c;
+
+        order_type m_q;
+
+        resizer_type m_resizer;
+
+    private:
+        void reset_iter(typename container_iterator_vector::iterator xiter_iter)
+        {
+            typename container_vector::iterator x_iter = m_xvec.begin();
+            while( x_iter != m_xvec.end() ) {
+                (*xiter_iter++) = (*x_iter++).begin();
+            }
+        }
+
+        void check_consitency()
+        {
+            const time_type* a_iter = &m_a[0];
+            const time_type* b_iter = &m_b[0];
+            const time_type* c_iter = &m_c[0];
+
+            unsigned short b_len = 1;
+            // check 1: a_i = sum_j b_ij 
+            while( a_iter != &m_a[+m_q-1] ) {
+                time_type tmp = 0.0;
+                const time_type* b_end = b_iter + b_len;
+                while( b_iter != b_end ) {
+                    tmp += *b_iter++;
+                }
+                b_len++;
+                if( *a_iter++ != tmp ) 
+                    throw butcher_tableau_consistency_exception();
+            }
+
+            // check 2: sum_i c_i * (a_i)^(k-1) = 1/k   k = 1..q
+            for( unsigned short k = 1; k <= m_q; k++ ) {
+                time_type tmp = 0.0;
+                a_iter = &m_a[0];
+                c_iter = &m_c[0];
+                if( k == 1 ) // special treatment for 0^0 = 1
+                    tmp += *c_iter++;
+                else
+                    c_iter++;
+                while( a_iter != &m_a[m_q-1] ) {
+                    tmp += (*c_iter++) * pow( *a_iter++ , k-1 );
+                }
+                if( std::abs(time_type(tmp - time_type(1.0)/time_type(k))) > 
+                    std::numeric_limits<time_type>::epsilon() ) {
+                    //std::clog << std::abs(time_type(tmp - time_type(1.0)/time_type(k))) << std::endl;
+                    throw butcher_tableau_order_condition_exception();
+                }
+            }
+        }
+
+
+    public:
+
+        /* Constructor
+           
+           Same as for the stepper_rk_generic class, but with a, b and c provided 
+           as time_type*. The order q of the integration scheme is given separately 
+           as parameter. a is a pointer to an array of time_type[] of length q-1,
+           b is of length 1 + 2 + ... q-1 = (q-1)(q-2)/2 and ordered as follows:
+           b[0] = b_21, b[1] = b_31, b[2] = b_32, b[3] = b_41, b[4] = b42 ...
+           c has length q.
+           
+        */
+        stepper_rk_generic_ptr( const time_type* a,
+                                const time_type* b,
+                                const time_type* c,
+                                const unsigned short q)
+            : m_a(a), m_b(b), m_c(c), m_q(q)
+        {
+            m_xvec.resize(m_q);
+            m_xiter_vec.resize(m_q);
+
+            check_consitency();
+        }
+
+        order_type order() const { return m_q; }
+
+        template< class DynamicalSystem >
+        void next_step( DynamicalSystem &system ,
+                        container_type &x ,
+                        container_type &dxdt ,
+                        time_type t ,
+                        time_type dt )
+        {
+            using namespace detail::it_algebra;
+            typename container_vector::iterator x_iter = m_xvec.begin();
+            typename container_iterator_vector::iterator xiter_iter = m_xiter_vec.begin();
+
+            (*x_iter) = dxdt;
+            (*xiter_iter++) = (*x_iter++).begin();
+
+            while( x_iter != m_xvec.end() ) {
+                m_resizer.adjust_size(x, (*x_iter));
+                (*xiter_iter++) = (*x_iter++).begin();
+            }
+            m_resizer.adjust_size(x, m_xtmp);
+            
+            x_iter = m_xvec.begin()+1;
+            
+            const time_type* a_iter = &m_a[0];
+            const time_type* b_iter = &m_b[0];
+            unsigned short b_len= 1;
+            while( x_iter != m_xvec.end() ) {
+                reset_iter(m_xiter_vec.begin());
+                const time_type* b_end = b_iter + b_len;
+                scale_sum_generic( m_xtmp.begin(), m_xtmp.end(),
+                                   b_iter, b_end, dt,
+                                   x.begin(), m_xiter_vec.begin() );
+                system( m_xtmp, *x_iter++ , t + dt*(*a_iter++) );
+                b_iter = b_end;
+                b_len++;
+            }
+
+            reset_iter(m_xiter_vec.begin());
+            scale_sum_generic( x.begin(), x.end(),
+                               &m_c[0], &m_c[m_q], dt,
+                               x.begin(), m_xiter_vec.begin() );
+        }
+
+
+
+        template< class DynamicalSystem >
+        void next_step( DynamicalSystem &system ,
+                        container_type &x ,
+                        time_type t ,
+                        time_type dt )
+        {
+            m_resizer.adjust_size(x, m_xtmp);
+            system(x, m_xtmp, t);
+            next_step( system, x, m_xtmp, t, dt);
+        }
+
+    };
+
 }
 }
 }

Modified: sandbox/odeint/libs/numeric/odeint/examples/lorenz_cmp_rk4_rk_generic.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/lorenz_cmp_rk4_rk_generic.cpp	(original)
+++ sandbox/odeint/libs/numeric/odeint/examples/lorenz_cmp_rk4_rk_generic.cpp	2009-11-16 17:18:44 EST (Mon, 16 Nov 2009)
@@ -21,7 +21,6 @@
 #include <tr1/array>
 
 #include <boost/numeric/odeint.hpp>
- // #include <boost/numeric/odeint/stepper_rk_generic.hpp>
 
 #define tab "\t"
 
@@ -76,13 +75,13 @@
     c[2] = 1.0/3.0;
     c[3] = 1.0/6.0;
 
-/*    const double a2[3] = {0.5, 0.5, 1.0 };
+    const double a2[3] = {0.5, 0.5, 1.0 };
     const double b2[6] = {0.5, 0.0, 0.5, 0.0, 0.0, 1.0};
-    const double c2[4] = {1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0};*/
+    const double c2[4] = {1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0};
 
     stepper_rk_generic< state_type > stepper_generic4(a, b, c);
 
-//    stepper_rk_generic_test< state_type > stepper_generic4_test(a2, b2, c2, 4);
+    stepper_rk_generic_ptr< state_type > stepper_generic4_ptr(a2, b2, c2, 4);
 
     clock_t start , end;
     double t;
@@ -117,6 +116,20 @@
     cout << "x after "<<olen<<" steps: "<<x2[0]<<tab<<x2[1]<<tab<<x2[2]<<endl;
     cout << "Time for "<<olen<<" steps: " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;
 
+    cout << endl << "###########################" << endl << endl;
+    cout << "Runge Kutta 4 via C-Array styled generic Runge Kutta implementation"<<endl;
+
+    t = 0.0;
+    start= clock();
+    for( size_t oi=1 ; oi<olen ; ++oi,t+=dt ) {
+        stepper_generic4_ptr.next_step( lorenz , x3 , t , dt );
+        if( oi < 5 )
+            cout << "x after step "<<oi<<":  "<<x3[0]<<tab<<x3[1]<<tab<<x3[2]<<endl;        
+    }
+    end = clock();
+    cout << "x after "<<olen<<" steps: "<<x3[0]<<tab<<x3[1]<<tab<<x3[2]<<endl;
+    cout << "Time for "<<olen<<" steps: " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;
+
 
     cout << endl << "###########################" << endl << endl;