[Boost-commit] svn:boost r72264 - sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/state_wrapper

Subject: [Boost-commit] svn:boost r72264 - sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/state_wrapper
From: mario.mulansky_at_[hidden]
Date: 2011-05-29 13:54:29

Author: mariomulansky
Date: 2011-05-29 13:54:29 EDT (Sun, 29 May 2011)
New Revision: 72264
URL: http://svn.boost.org/trac/boost/changeset/72264

Log:
crude state wrapper implementation
Added:
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/state_wrapper/
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/state_wrapper/Jamfile   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/state_wrapper/explicit_euler.hpp   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/state_wrapper/test_gsl_vector.cpp   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/state_wrapper/test_vector.cpp   (contents, props changed)

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/state_wrapper/Jamfile
==============================================================================
Binary file. No diff available.

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/state_wrapper/explicit_euler.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/state_wrapper/explicit_euler.hpp	2011-05-29 13:54:29 EDT (Sun, 29 May 2011)
@@ -0,0 +1,72 @@
+#include <iostream>
+
+#include <boost/range.hpp>
+
+#include <boost/numeric/odeint/algebra/detail/for_each.hpp>
+#include <boost/numeric/odeint/algebra/default_operations.hpp>
+
+template< class V >
+struct state_wrapper
+{
+    typedef typename V::value_type value_type;
+
+    V &m_v;
+    V m_v_;
+
+    state_wrapper( V &v ) : m_v( v ) { }
+
+    state_wrapper() : m_v( m_v_ ) { }
+
+    typename boost::range_iterator<V>::type begin()
+    { return boost::begin( m_v ); }
+
+    typename boost::range_iterator<V>::type end()
+    { return boost::end( m_v ); }
+
+    void resize( V &v )
+    {
+        m_v.resize( boost::size( v ) );
+    }
+
+    bool same_size( V &v )
+    { return ( boost::size( m_v ) == boost::size( v ) ); }
+};
+
+
+struct range_algebra
+{
+    template< class S1 , class S2 , class S3 , class Op >
+    static void for_each3( S1 &s1 , S2 &s2 , S3 &s3 , Op op )
+    {
+        boost::numeric::odeint::detail::for_each3( s1.begin() , s1.end() , s2.begin() , s3.begin() , op );
+    }
+};
+
+template< typename StateType >
+class explicit_euler {
+
+public:
+
+    typedef double value_type;
+    typedef double time_type;
+    typedef StateType state_type;
+    typedef state_wrapper< state_type > wrapped_state_type;
+
+    explicit_euler() { };
+
+    template< class System , class StateInOut >
+    void do_step( System system , StateInOut &inout , const time_type &t , const time_type &dt )
+    {
+        state_wrapper< StateInOut > x( inout );
+
+        if( !m_dxdt.same_size( inout ) )
+            m_dxdt.resize( inout );
+
+        system( inout , m_dxdt.m_v , t );
+
+        range_algebra::for_each3( x , x , m_dxdt , typename boost::numeric::odeint::default_operations::template scale_sum2< value_type , time_type >( 1.0 , dt ) );
+    }
+
+private:
+    wrapped_state_type m_dxdt;
+};

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/state_wrapper/test_gsl_vector.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/state_wrapper/test_gsl_vector.cpp	2011-05-29 13:54:29 EDT (Sun, 29 May 2011)
@@ -0,0 +1,89 @@
+#include <iostream>
+#include <gsl/gsl_vector.h>
+
+#include "explicit_euler.hpp"
+
+using namespace std;
+
+typedef gsl_vector *state_type;
+
+const double sigma = 10.0;
+const double R = 28.0;
+const double b = 8.0 / 3.0;
+
+void lorenz( const state_type x , state_type dxdt , double t )
+{
+    gsl_vector_set( dxdt , 0 , sigma * ( gsl_vector_get(x , 1 ) - gsl_vector_get( x , 0 ) ) );
+    gsl_vector_set( dxdt , 1 , R * gsl_vector_get( x , 0 ) - gsl_vector_get( x , 1 ) - gsl_vector_get( x , 0 ) * gsl_vector_get( x , 2) );
+    gsl_vector_set( dxdt , 2 , gsl_vector_get( x , 0 ) * gsl_vector_get( x , 1 ) - b * gsl_vector_get( x , 2) );
+}
+
+template<>
+struct state_wrapper< state_type >
+{
+    typedef double value_type;
+
+    state_type m_v;
+
+    state_wrapper( state_type &v ) : m_v( v ) {
+        cout << m_v->size << endl;
+    }
+
+    state_wrapper( )
+    {
+        m_v->owner = 0;
+        m_v->size = 0;
+        m_v->stride = 0;
+        m_v->data = 0;
+        m_v->block = 0;
+    }
+
+    double* begin()
+    { return m_v->data; }
+
+    double* end()
+    { return m_v->data + m_v->size; }
+
+    void resize( state_type &v )
+    {
+        cout << v->size << " " << m_v->owner << " " << v->owner << endl;
+
+        if( m_v->owner != 0 )
+        {
+            gsl_block_free( m_v->block );
+        }
+        m_v->size = 0;
+
+        cout << v->size << endl;
+
+        if( v->size == 0 ) return;
+
+        gsl_block *block = gsl_block_alloc( v->size );
+        if( block == 0 ) throw std::bad_alloc( );
+
+        m_v->data = block->data ;
+        m_v->size = v->size;
+        m_v->stride = 1;
+        m_v->block = block;
+        m_v->owner = 1;
+    }
+
+    bool same_size( state_type &v )
+    { return ( m_v->size == v->size ); }
+};
+
+int main() {
+
+    explicit_euler< state_type > euler;
+
+    state_type x = gsl_vector_alloc( 3 );
+    gsl_vector_set( x , 0 , 1.0);
+    gsl_vector_set( x , 1 , 1.0);
+    gsl_vector_set( x , 2 , 2.0);
+
+    euler.do_step( lorenz , x , 0.0 , 0.1 );
+
+    cout << gsl_vector_get( x , 0 ) << "  " << gsl_vector_get( x , 1 ) << "  " << gsl_vector_get( x , 2 ) << endl;
+
+    gsl_vector_free( x );
+}

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/state_wrapper/test_vector.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/state_wrapper/test_vector.cpp	2011-05-29 13:54:29 EDT (Sun, 29 May 2011)
@@ -0,0 +1,31 @@
+#include <iostream>
+#include <vector>
+
+#include "explicit_euler.hpp"
+
+using namespace std;
+
+typedef vector< double > state_type;
+
+const double sigma = 10.0;
+const double R = 28.0;
+const double b = 8.0 / 3.0;
+
+void lorenz( const state_type &x , state_type &dxdt , 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];
+}
+
+int main()
+{
+    explicit_euler< state_type > euler;
+
+    state_type x(3);
+    x[0] = 1.0; x[1] = 1.0; x[2] = 2.0;
+
+    euler.do_step( lorenz , x , 0.0 , 0.1 );
+
+    cout << x[0] << "  " << x[1] << "  " << x[2] << endl;
+}