
Comments and Discussions



Hi, thanks for providing this library.
I am simulating an electrical circuit using odeint. The results are not compatible with matlab results (using ode23). I have tried fixed step size adaptive stpe size and dense output. I have tried reducing the step sizes, and error tolerances. I know that matlab results are correct because I simulated my circuit in a dedicated circuit simulator, and matlab matches that.
Would you have any suggestions on debugging the problem? If needed, I can provide my code, which is not large.
Thanks,





Hi,
it would be nice to have an example on how to integrate nonautonomous systems via thrust. The problem is how to pass the time in the functor.
I can't see a simple way to do it.
Can somebody help me?
Thanks for this library.
Michele





Hi, this is pretty easy. Just pass the time to the functor you pass to thrusts for_each:
struct ode_functor
{
double m_t;
ode_functor( double t ) : m_t( t ) { }
template< class T >
__host__ __device__
void operator()( T t ) const
{
thrust::get< 1 > =  thrust::get< 0 > * m_t;
}
};
struct ode
{
template< class State , class Deriv >
void operator()( const State &x , Deriv &dxdt , value_type t ) const
{
thrust::for_each(
thrust::make_zip_iterator( thrust::make_tuple(
boost::begin( x ) ,
boost::begin( dxdt ) ) ),
thrust::make_zip_iterator( thrust::make_tuple(
boost::end( x )
boost::end( dxdt ) ) ),
ode_functor( t ) );
}
};





Many many thanks. Exactly what I needed. I'm not very expert about c++ so I didn't understand how really Thrust works with Odeint. Now I did a step forward.
How about including an example like this in the odeint tutorial? It's already a very good work but given that all the examples are about autonomous systems, it can be more effective with a simple example like this.
Thanks,
Michele
modified 17Mar13 12:06pm.





Thanks for that hint. We will try to include nonautonomous system in the Thrust tutorial.





Hello,
great library, respect
I have a problem with negative timesteps. All functions are designed for positive timesteps. Can you tell me, what I' must changing?
Thanks, Thomas






Thank you. The actuall version works well





Hi,
sorry to disturb but I am new to the code project and excited about using odeint v2 for academic research. I am trying to use the implicit euler stepper for the first step in a system of nonlinear ODEs and when I type stepper.do_step (system, x, r , dr) I get an error message saying the the function do_step can not be found. I thought that this was a function every stepper had. Am I doing something wrong? I tried looking for an example with the implicit euler used but I have not seen any.
Any help would be appreciated. Thank you in advance!





Hi,
maybe you are misusing the stepper method of the implict_euler. As system function it expects a pair of function object, one returning the r.h.s. of the ode and the other one the Jacobian of the r.h.s. An example can be found here (ok, it is the rosenbrock stepper, but this one is very similar to the implicit euler).
Hope this helps,
Karsten





Hi,
thank you for your quick reply. I finally got to try the example you gave me with implicit euler but I still have problems. To be honest, I'm pretty new to template usage and probably my mistakes are trivial but I still need help. I just need to do the first step in my problem implicitly and use runge kutta for the rest. Below I have written the code I have so far, using the stiff system from your web page. Please, let me know what I am doing wrong!
#include <iostream>
#include <fstream>
#include <utility>
#include <boost/numeric/odeint.hpp>
#include <boost/spirit/include/phoenix.hpp>
typedef std::vector< double > state_type;
using namespace std;
using namespace boost::numeric::odeint;
namespace phoenix = boost::phoenix;
//[ stiff_system_definition
typedef boost::numeric::ublas::vector< double > vector_type;
typedef boost::numeric::ublas::matrix< double > matrix_type;
struct stiff_system
{
void operator()( const vector_type &x , vector_type &dxdt , double /* t */ )
{
dxdt[ 0 ] = 101.0 * x[ 0 ]  100.0 * x[ 1 ];
dxdt[ 1 ] = x[ 0 ];
}
};
struct stiff_system_jacobi
{
void operator()( const vector_type & /* x */ , matrix_type &J , const double & /* t */ , vector_type &dfdt )
{
J( 0 , 0 ) = 101.0;
J( 0 , 1 ) = 100.0;
J( 1 , 0 ) = 1.0;
J( 1 , 1 ) = 0.0;
dfdt[0] = 0.0;
dfdt[1] = 0.0;
}
};
int main( int argc , char **argv )
{
vector_type x( 3 , 1.0 );
implicit_euler< double, vector_type > doit;
doit.do_step( make_pair( stiff_system() , stiff_system_jacobi() ) , x , 0.0, 0.05);
system("pause");
return 0;
}





Hi,
I looked at your example and found two little points:
1. the implicit_euler is instantiated only with one template argument double
2. the Jacobian is calculated without explicitly calculating the partial derivative of the system function with respect to the time. Simply replace stiff_system_jacobi with
struct stiff_system_jacobi
{
void operator()(const vector_type &x , matrix_type &J , const double &t )
{
J(0, 0) = 101.0;
J(0, 1) = 100.0;
J(1, 0) = 1.0;
J(1, 1) = 0.0;
}
};
and everything should work as expected.
There are some small inconsistencies with the implicit steppers. The Rosenbrock stepper (another implicit stepper) expect the partial derivative of the system function w.r.t. the time explicitly, while the implicit Euler does not.





Hi,
thank you for your help. I opted to use the Rosenbrock method since my system of ODE has explicit t dependence. However, I still have problems. This time it seems to be due to division by t^2 at t=0 which is the beginning of my domain. I thought that since it is an implicit method it should deal wtih this kind of situations. The message I get is "This application has requested the Runtime to terminate it in an unusual way. Please Contact the application's support team for more information." And my code is below:
#include <iostream>
#include <iomanip>
#include <fstream>
#include <utility>
#include <boost/numeric/odeint.hpp>
#include <boost/spirit/include/phoenix.hpp>
using namespace std;
using namespace boost::numeric::odeint;
namespace phoenix = boost::phoenix;
typedef boost::numeric::ublas::vector< double > vector_type;
typedef boost::numeric::ublas::matrix< double > matrix_type;
vector_type x( 2 );
vector_type dxdt ( 2 );
void stiffsystem( vector_type &x , vector_type &dxdt , double t )
{
dxdt[ 0 ] = 101.0 * x[ 0 ] * x[ 1 ]/ pow(t,2.0);
dxdt[ 1 ] = x[ 0 ]* pow(t,2.0);
}
struct stiff_system
{
void operator()( const vector_type &x , vector_type &dxdt , double t )
{
dxdt[ 0 ] = 101.0 * x[ 0 ] * x[ 1 ]/pow(t,2.0);
dxdt[ 1 ] = x[ 0 ]*pow(t,2.0);
}
};
struct stiff_system_jacobi
{
void operator()( const vector_type & x , matrix_type &J , const double & t , vector_type &dfdt )
{
J( 0 , 0 ) = 101.0* x[1]/pow(t,2.0);
J( 0 , 1 ) = 100.0*x[0]/pow(t,2.0);
J( 1 , 0 ) = pow(t,2.0);
J( 1 , 1 ) = 0.0;
dfdt[0] = 2.0*101.0 * x[ 0 ] * x[ 1 ]*pow(t,3.0);
dfdt[1] = 2.0*t*x[0];
}
};
struct streaming_observer
{
std::ofstream &m_out;
streaming_observer( std::ofstream &out) : m_out( out ) {}
void operator()( const vector_type &x, double t ) const
{
m_out << t;
for( size_t i=0 ; i < x.size() ; ++i )
{
m_out << "\t" << x[i];
}
m_out << "\n";
}
};
int main( int argc , char **argv )
{
ofstream output;
ifstream input;
output.open("data.txt");
const double dt = 0.01;
double t;
double inv[] = {1.0 ,0.5};
dxdt [0] = 0.0;
dxdt [1] = 0.0;
x[0] = inv[0];
x[1] = inv[1];
integrate_adaptive( make_dense_output< rosenbrock4< double > >( 1.0e6 , 1.0e6 ) ,
make_pair( stiff_system() , stiff_system_jacobi() ) ,
x , 0.0 , 0.1 , dt ,
streaming_observer(output)),
cout << x[0] << " " << x[1] <<" "<<endl;
output.close();
input.open ("data.txt");
while(!input.eof())
{
input >> t >> x[0] >> x[1];
};
input.close();
cout << t << " " << x[0] << " "<< x[1] <<" "<<endl;
output.open ("data.txt", fstream::app);
runge_kutta4< vector_type > rk4;
for( size_t i=0 ; i<1000 ; i++)
{
t+=dt;
rk4.do_step( stiffsystem , x , t , dt );
if (abs(x[0]) < 0.1)
break;
output << t << " " << x[0] <<"\t"<< x[1] << "\n";
}
output.close();
cout << t << " " << x[0] << " "<< x[1] <<" "<<endl;
system("pause");
return 0;
}
Please, help! I really need this bit to be working as it is the main reason behind using the code you provide here instead of doing the first steps by hand. Again, thank you!





Hi, in your system you are dividing by t^2. And your starting time is 0, such that a term 1 / 0 ^2 appears in your system definition. If I change the start time from 0.0 to 0.01 everything seems fine.
integrate_adaptive( make_dense_output< rosenbrock4< double > >( 1.0e6 , 1.0e6 ) ,
make_pair( stiff_system() , stiff_system_jacobi() ) ,
x , 0.01 , 0.1 , dt ,
streaming_observer(output) );
Btw. you can use t*t for pow(t,2.0) which should be much faster.





Hi,
I think I already mentioned that I'm starting at 0. But that should be a problem only for explicit methods whcih use the current point to go forward. Impicit methods like backward Euler should use the t+dt point to find the value of the function at t+dt, which is the main reason why I want to use your code. I thought Rosenbrock had the same capbility. Am I wrong? Is there a way to do that? I have already tried using 0.01 instead of 0 but that leads to big numerical errors later in calculations. Thank you!





Ok, Rosenbrock is a semiimplicit method. It calculates the Jacobian at t=0 and evaluates the function at t=0. I think you can use implicit euler method. But it expects the system function in a slightly different way, such that you have to rewrite it. But this should only a minor change for your problem.
modified 21Jun12 3:46am.





Hi,
I did the changes you have previously mentioned but neither do step nor integrate_adaptive seem to work with the Euler method. You mentioned that I have to rewrite the function. Would you mind showing me how? There's no example on how to use implicit_euler. My changed code is below:
#include <iostream>
#include <iomanip>
#include <fstream>
#include <utility>
#include <boost/numeric/odeint.hpp>
#include <boost/spirit/include/phoenix.hpp>
using namespace std;
using namespace boost::numeric::odeint;
namespace phoenix = boost::phoenix;
// stiff_system_definition
typedef boost::numeric::ublas::vector< double > vector_type;
typedef boost::numeric::ublas::matrix< double > matrix_type;
typedef vector< double > state_type;
const double PI = 4.134;
const double k = 1.0;
const double lambda = 1.0;
const double gamma = 2.0;
const double g = 1.0;
state_type x( 2 );
double Lane_Emden (double rho, double a, double b)
{
return (a*pow(rho,b));
}
double d_Lane_Emden (double rho, double a, double b)
{
return (a*b*pow(rho,b1.0));
}
double dd_Lane_Emden (double rho, double a, double b)
{
return (a*b*(b1.0)*pow(rho,b2.0));
}
void stiffsystemtest( state_type &x , state_type &dxdt , double t )
{
dxdt[ 0 ] =  g * x[ 0 ] * x[ 1 ] / pow(t,2.0) / d_Lane_Emden(x[0],k,gamma);
dxdt[ 1 ] = 4 * PI * pow(t,2.0)* x[ 0 ];
};
struct stiff_systemtest
{
template< class State >
void operator()( const State &x , State &dxdt , double t )
{
dxdt[ 0 ] =  g * x[ 0 ] * x[ 1 ] / pow(t,2.0) / d_Lane_Emden(x[0],k,gamma);
dxdt[ 1 ] = 4.0 * PI * pow(t,2.0)* x[ 0 ];
}
};
struct stiff_system_jacobitest
{
template< class State , class Matrix >
void operator()( const State &x , Matrix &J , const double &t)
{
J( 0 , 0 ) =  g * x[ 1 ] / pow(t,2.0) / d_Lane_Emden(x[0],k,gamma)
+ g * x[ 0 ] * x[ 1 ] / pow(t,2.0)
/ pow(d_Lane_Emden(x[0],k,gamma),2.0) * dd_Lane_Emden(x[0],k,gamma);
J( 0 , 1 ) =  g * x[ 0 ] / pow(t,2.0) / d_Lane_Emden(x[0],k,gamma);
J( 1 , 0 ) = 4 * PI * pow(t,2.0);
J( 1 , 1 ) = 0.0;
//dfdt[0] = 0.0;
//dfdt[1] = 0.0;
}
};
struct streaming_observer
{
std::ofstream &m_out;
streaming_observer( std::ofstream &out) : m_out( out ) {}
void operator()( const vector_type &x, double t ) const
{
m_out << t;
for( size_t i=0 ; i < x.size() ; ++i )
{
m_out << "\t" << x[i];
}
m_out << "\n";
}
};
int main( int argc , char **argv )
{
ofstream output;
ifstream input;
output.open("testdata.txt");
const double dt = 0.01;
double t;
double inv[] = {lambda ,0.0};
x[0] = inv[0];
x[1] = inv[1];
//implicit_euler< double > doit;
//doit.do_step( make_pair( stiff_systemtest() , stiff_system_jacobitest() ) , x , 0.0, 0.05);
integrate_adaptive(implicit_euler< double >( 1.0e6 ) ,
make_pair( stiff_systemtest() , stiff_system_jacobitest() ) ,
x , 0.0 , 0.01 , 0.1, streaming_observer(output));
output.close();
input.open ("testdata.txt");
while(!input.eof())
{
input >





The trick is to put the time into the definition of the system. So your system consists of three equation, where one describes trivially the time:
struct stiff_system
{
void operator()( const vector_type &x , vector_type &dxdt , double t )
{
dxdt[ 0 ] = 101.0 * x[ 0 ] * x[ 1 ] / t / t;
dxdt[ 1 ] = x[ 0 ] * t * t;
dxdt[ 2 ] = 1.0;
}
};
struct stiff_system_jacobi
{
void operator()( const vector_type & x , matrix_type &J , const double & t )
{
J( 0 , 0 ) = 101.0 * x[1] / t / t;
J( 0 , 1 ) = 101.0 * x[0] / t / t;
J( 0 , 2 ) = 202.0 * x[0] * x[1] / t / t / t;
J( 1 , 0 ) = t * t ;
J( 1 , 1 ) = 0.0;
J( 1 , 2 ) = 2 * t * x[0];
J( 2 , 0 ) = 0.0;
J( 2 , 1 ) = 0.0;
J( 2 , 2 ) = 0.0;
}
};
Then you can use it, like this
const double dt = 0.001;
vector_type x( 3 );
x[0] = 1.0;
x[1] = 0.5;
x[2] = 0.0;
integrate_adaptive( implicit_euler< double >() ,
make_pair( stiff_system() , stiff_system_jacobi() ) ,
x , 0.0 , 0.1 , dt ,
streaming_observer( cout ) );
Do the result fullfil your expectation?
Note, that the implicit_euler does not have any dense output or step size control, such that integrate_adaptive is only an integrate_const. Sorry for the inconsistency of the system functions interface. Unifying the implicit steppers system function is an important point on our ToDo list. Of course, we also plan to implement more implicit methods, but it will take some time.





Hi,
I used the implicit euler stepper together with runge kutta 4 and they worked like a charm! Thank you for the help provided. By the way, I think it would be a good idea to add examples on how each stepper works since there some differences in their implementations. For instance, integrateadaptive shouldn't be used at all with implicit Euler since there is integrate_const which gives better results. I might bug you again in the near future!





Ok, thank you for the hint. We will add an example about the implicit_euler in the tutorial section.
By the way, in the case of the implicit_euler integrate_adaptive and integrate_const should give the same results. Of course, integrate_adaptive is misleading since it simply redirects in this case to integrate_adaptive.





Hello,
I had problems when I tried to start the example file Lorenz_gmpxx.cpp. I am a freshman for programming but I have to solve ODEs in my research. Here is what I did,
1. I download boost_1_49_0, gmp_5.0.4 and installed them in this directory /usr/local/
2. export BOOST_BOOT=/use/local/boost_1_49_0
3. I tried to run bjam under the directory odeintV2, showed
...failed updating 6 targets...
...skipped 122 targets...
then I run sudo bjam, showed
Unable to load Boost.Build: could not find "boost build.jam"

Attempted search from /home/mingzi/odeintv2 up to the root
at /usr/local/share/boostbuild
and in these directories from BOOST_BUILD_PATH and BOOST_ROOT: /usr/share/boostbuild.
Please consult the documentation at 'http:
4. I tried to run
g++ I /usr/local/boost_1_49_0 I /usr/local/gmp5.0.4/ I /home/mingzi/odeintv2/ lorenz_gmpxx.cpp o lorenz
and it shows(just a part of the message):
collect2: ld returned 1 exit status
<pre lang="vb">/tmp/ccAnjsOr.o: In function `main':
lorenz_gmpxx.cpp text+0x26): undefined reference to `__gmpf_set_default_prec'
/tmp/ccAnjsOr.o: In function `__gmp_unary_minus::eval(__mpf_struct*, __mpf_struct const*)':
lorenz_gmpxx.cpp text._ZN17__gmp_unary_minus4evalEP12__mpf_structPKS0_[__gmp_unary_minus::eval(__mpf_struct*, __mpf_struct const*)]+0x14): undefined reference to `__gmpf_neg'
/tmp/ccAnjsOr.o: In function `__gmp_binary_plus::eval(__mpf_struct*, __mpf_struct const*, __mpf_struct const*)':
lorenz_gmpxx.cpp text._ZN17__gmp_binary_plus4evalEP12__mpf_structPKS0_S3_[__gmp_binary_plus::eval(__mpf_struct*, __mpf_struct const*, __mpf_struct const*)]+0x1b): undefined reference to `__gmpf_add'</pre>
Can you please help me to figure out where the problem is? and any suggestions will be appreciated.
Greetings,
Mingzi





Stupid question: in your is a typo:
export BOOST_BOOT=/use/local/boost_1_49_0
is this typo in your .bashrc? It would explain the errors in point 3.
For point 4 you need to link against gmpxx and gmp, simply add
lgmpxx lgmp
to your command:
g++ I /usr/local/boost_1_49_0 I /usr/local/gmp5.0.4/ I /home/mingzi/odeintv2/ lorenz_gmpxx.cpp o lorenz lgmpxx lgmp





Hello
yesterday i downloaded the odeintv2 package, but i have no clue how to compile or setup it.
I have a pc with linux (opensuse 12.1) and boost 1.461.
I know the classic unixtools like make, gcc and friends, (also c and some c++) but i do not know how to use bjam,
bootstrap and this. i googled around, reading the boost site but this did not help me.
Can one friendly soul give me a hint how to start to setup odeintv2 from zero?
greetings Michael





There is nothing to setup. odeint is completely header only, all you have to do is to specify the directory during compilation. For gcc this is easy, just add the compiler flags I/path/to/odeint I/path/to_boost.
If you want to compile the examples and unit tests of odeint you have to
* download boost
* create an environment variable BOOST_ROOT pointing to the boost directory
* run bjam in the odeint directory
If bjam is not available you can
* run booststrap.sh n BOOST_ROOT (which creates bjam here)
* copy bjam in a directory like /usr/local/bin or $HOME/bin such that it can be found anywhere on your system
HTH





thanks for the answer. What i did:
export BOOST_ROOT=/usr/share/doc/packages/boost1.46.1 ((Here boost is installed when using the opensuse dvd))
# ~/downloads/odeintv2> bjam
Unable to load Boost.Build: could not find "boostbuild.jam"

Attempted search from /home/mma/downloads/odeintv2 up to the root
at /usr/share/boostbuild
and in these directories from BOOST_BUILD_PATH and BOOST_ROOT: /usr/share/boostbuild, /usr/share/doc/packages/boost1.46.1.
Please consult the documentation at 'http://www.boost.org'.
i fear there is lacking something in the installation. find only finds:
# find /usr/share/doc/packages/boost1.46.1 name boostbuild.jam
/usr/share/doc/packages/boost1.46.1/libs/python/example/quickstart/boostbuild.jam
no other boostbuild.jam. Is this possible ?
regards Michael





Yes, this possible. The installation from opensuse does not install the complete boost tree. I would suggest that you download boost and place it into you home directory.





Hello
i succeeded to setup odeintv2 and could start all example
But when i tried to compile the thrust examples for using cuda it failed.
(I have installed cuda 4.1, gcc3.6, thrust 1.48.0 and 2 nvidia gtx580 cards):
mma@arnie /downloads/odeintv2/libs/numeric/odeint/examples/thrust> make
/usr/local/cuda/bin/nvcc I/usr/share/doc/packages/boost_1_48_0 I I/usr/local/cuda/include I../../../../.. arch sm_20 compilerbindir=/usr/bin/g++ Xcompiler fopenmp DTHRUST_DEVICE_BACKEND=THRUST_DEVICE_BACKEND_OMP o phase_oscillator_chain.co c phase_oscillator_chain.cu
In file included from ../../../../../boost/numeric/odeint.hpp:21 ,
from phase_oscillator_chain.cu:31:
../../../../../boost/numeric/odeint/config.hpp:17 : warning: "FUSION_MAX_VECTOR_SIZE" redefined [enabled by default]
and so on and finally
instantiation of "size_t boost::numeric::odeint::integrate_const(Stepper, System, State &, Time, Time, Time) [with Stepper=boost::numeric::odeint::runge_kutta4, System=phase_oscillators, State=state_type, Time=double]"
phase_oscillator_chain.cu(156): here
25 errors detected in the compilation of "/tmp/tmpxft_00006715_000000004_phase_oscillator_chain.cpp1.ii".
make: *** [phase_oscillator_chain.co] Fehler 2
any ideas?
greetings Michael





I am sorry, I have no idea. Maybe it does not work with sm_20, maybe you need a more recent version of gcc? Can you try arch sm_13? Can you show some specific error messages?





Hello
here is the check_thrust.cu example.
I remember the Nvidia SDK said to use gcc4.4 but i managed to compile all nvidias examples with gcc3.6 and all examples i tested were working
~/downloads/odeintv2/libs/numeric/odeint/test_external/thrust>
make
/usr/local/cuda/bin/nvcc O3 I/usr/share/doc/packages/boost_1_48_0 I I/usr/local/cuda/include I../../../../.. compilerbindir=/usr/bin/g++ o check_thrust.co c check_thrust.cu
../../../../../boost/numeric/odeint/external/thrust/thrust_resize.hpp(46): error: a template argument list is not allowed in a declaration of a primary template
../../../../../boost/numeric/odeint/external/thrust/thrust_resize.hpp(55): error: a template argument list is not allowed in a declaration of a primary template
../../../../../boost/numeric/odeint/external/thrust/thrust_resize.hpp(64): error: a template argument list is not allowed in a declaration of a primary template
/usr/share/doc/packages/boost_1_48_0/boost/range/size.hpp(32): warning: controlling expression is constant
detected during:
instantiation of "boost::range_difference::type boost::range_detail::range_calculate_size(const SinglePassRange &) [with SinglePassRange=state_type]"
(47): here
instantiation of "boost::range_difference::type boost::size(const SinglePassRange &) [with SinglePassRange=state_type]"
../../../../../boost/numeric/odeint/util/state_wrapper.hpp(50): here
instantiation of "__nv_bool boost::numeric::odeint::state_wrapper::same_size(const StateIn &) const [with V=state_type, StateIn=state_type]"
../../../../../boost/numeric/odeint/util/state_wrapper.hpp(57): here
instantiation of "__nv_bool boost::numeric::odeint::state_wrapper::resize(const StateIn &) [with V=state_type, StateIn=state_type]"
../../../../../boost/numeric/odeint/util/resizer.hpp(31): here
instantiation of "__nv_bool boost::numeric::odeint::adjust_size_by_resizeability(ResizeState &, const State &, boost::true_type) [with ResizeState=boost::numeric::odeint::state_wrapper, State=state_type]"
../../../../../boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp(146): here
instantiation of "__nv_bool boost::numeric::odeint::explicit_stepper_base::resize(const StateIn &) [with Stepper=boost::numeric::odeint::euler, Order=(unsigned short)1U, State=state_type, Value=base_type, Deriv=state_type, Time=base_type, Algebra=boost::numeric::odeint::thrust_algebra, Operations=boost::numeric::odeint::thrust_operations, Resizer=boost::numeric::odeint::initially_resizer, StateIn=state_type]"
../../../../../boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp(173): here
instantiation of "void boost::numeric::odeint::explicit_stepper_base::do_step_v1(System, StateInOut &, const Time &, const Time &) [with Stepper=boost::numeric::odeint::euler, Order=(unsigned short)1U, State=state_type, Value=base_type, Deriv=state_type, Time=base_type, Algebra=boost::numeric::odeint::thrust_algebra, Operations=boost::numeric::odeint::thrust_operations, Resizer=boost::numeric::odeint::initially_resizer, System=void (*)(const state_type &, state_type &, base_type), StateInOut=state_type]"
../../../../../boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp(95): here
instantiation of "void boost::numeric::odeint::explicit_stepper_base::do_step(System, StateInOut &, const Time &, const Time &) [with Stepper=boost::numeric::odeint::euler, Order=(unsigned short)1U, State=state_type, Value=base_type, Deriv=state_type, Time=base_type, Algebra=boost::numeric::odeint::thrust_algebra, Operations=boost::numeric::odeint::thrust_operations, Resizer=boost::numeric::odeint::initially_resizer, System=void (*)(const state_type &, state_type &, base_type), StateInOut=state_type]"
check_thrust.cu(45): here
instantiation of "void check_stepper_concept(Stepper &, System, Stepper::state_type &) [with Stepper=boost::numeric::odeint::euler, System=void (*)(const state_type &, state_type &, base_type)]"
check_thrust.cu(58): here
3 errors detected in the compilation of "/tmp/tmpxft_00007105_000000004_check_thrust.cpp1.ii".
greetings Michael





Hmm, I have no idea. Can you try with gcc 4.4?





Hello
i tried with gcc4.3. > The same errors
But now i suspect "thrust or boost its too new". Wich thrust and boost versions did you use?
Greetings Michael





I use gcc 4.4, boost 1.45 and thrust 1.4(?). I think the concrete boost or thrust version is not so important. Maybe gcc 4.4 would solve the problem.





I got the same report with vs2008, odeint v2, boost 1.45 and thrust 1.4 as follows,
1> Build started: Project: template, Configuration: Debug Win32 
1>Compiling with CUDA Build Rule...
1>"H:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v4.0\\bin\nvcc.exe" gencode=arch=compute_10,code=\"sm_10,compute_10\" gencode=arch=compute_20,code=\"sm_20,compute_20\" machine 32 ccbin "H:\Program Files (x86)\Microsoft Visual Studio 9.0\VC\bin" Xcompiler "/EHsc /W3 /nologo /O2 /Zi /MT " I"H:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v4.0\/include" I"./" I"../../common/inc" I"../../../shared/inc" I"H:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v4.0\\include" maxrregcount=0 compile o "Debug/main.cu.obj" main.cu
1>main.cu
1>E:\odeintv2\boost/numeric/odeint/external/thrust/thrust_resize.hpp(46): error: a template argument list is not allowed in a declaration of a primary template
1>E:\odeintv2\boost/numeric/odeint/external/thrust/thrust_resize.hpp(55): error: a template argument list is not allowed in a declaration of a primary template
1>E:\odeintv2\boost/numeric/odeint/external/thrust/thrust_resize.hpp(64): error: a template argument list is not allowed in a declaration of a primary template
1>E:\odeintv2\boost/range/size.hpp(29): warning: controlling expression is constant
1> detected during:
1> instantiation of "boost::range_difference::type boost::size(const T &) [with T=state_type]"
1>E:\odeintv2\boost/numeric/odeint/util/state_wrapper.hpp(50): here
1> instantiation of "bool boost::numeric::odeint::state_wrapper::same_size(const StateIn &) const [with V=state_type, StateIn=state_type]"
1>E:\odeintv2\boost/numeric/odeint/util/state_wrapper.hpp(57): here
1> instantiation of "bool boost::numeric::odeint::state_wrapper::resize(const StateIn &) [with V=state_type, StateIn=state_type]"
1>E:\odeintv2\boost/numeric/odeint/util/resizer.hpp(31): here
1> instantiation of "bool boost::numeric::odeint::adjust_size_by_resizeability(ResizeState &, const State &, boost::true_type) [with ResizeState=boost::numeric::odeint::state_wrapper, State=state_type]"
1>E:\odeintv2\boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp(146): here
1> instantiation of "bool boost::numeric::odeint::explicit_stepper_base::resize(const StateIn &) [with Stepper=boost::numeric::odeint::euler, Order=(unsigned short)1U, State=state_type, Value=base_type, Deriv=state_type, Time=base_type, Algebra=boost::numeric::odeint::thrust_algebra, Operations=boost::numeric::odeint::thrust_operations, Resizer=boost::numeric::odeint::initially_resizer, StateIn=state_type]"
1>E:\odeintv2\boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp(173): here
1> instantiation of "void boost::numeric::odeint::explicit_stepper_base::do_step_v1(System, StateInOut &, const Time &, const Time &) [with Stepper=boost::numeric::odeint::euler, Order=(unsigned short)1U, State=state_type, Value=base_type, Deriv=state_type, Time=base_type, Algebra=boost::numeric::odeint::thrust_algebra, Operations=boost::numeric::odeint::thrust_operations, Resizer=boost::numeric::odeint::initially_resizer, System=void (*)(const state_type &, state_type &, base_type), StateInOut=state_type]"
1>E:\odeintv2\boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp(95): here
1> instantiation of "void boost::numeric::odeint::explicit_stepper_base::do_step(System, StateInOut &, const Time &, const Time &) [with Stepper=boost::numeric::odeint::euler, Order=(unsigned short)1U, State=state_type, Value=base_type, Deriv=state_type, Time=base_type, Algebra=boost::numeric::odeint::thrust_algebra, Operations=boost::numeric::odeint::thrust_operations, Resizer=boost::numeric::odeint::initially_resizer, System=void (*)(const state_type &, state_type &, base_type), StateInOut=state_type]"
1>h:/ProgramData/NVIDIA Corporation/NVIDIA GPU Computing SDK 4.0/C/src/template_my/main.cu(286): here
1> instantiation of "void check_stepper_concept(Stepper &, System, Stepper::state_type &) [with Stepper=boost::numeric::odeint::euler, System=void (*)(const state_type &, state_type &, base_type)]"
1>h:/ProgramData/NVIDIA Corporation/NVIDIA GPU Computing SDK 4.0/C/src/template_my/main.cu(299): here
1>3 errors detected in the compilation of "H:/Users/wang/AppData/Local/Temp/tmpxft_00010e20_000000008_main.compute_10.cpp1.ii".
1>Project : error PRJ0019: A tool returned an error code from "Compiling with CUDA Build Rule..."
1>Build log was saved at "file://h:\ProgramData\NVIDIA Corporation\NVIDIA GPU Computing SDK 4.0\C\src\template_my\Debug\BuildLog.htm"
1>template  4 error(s), 1 warning(s)
========== Build: 0 succeeded, 1 failed, 1 uptodate, 0 skipped ==========





Hello
in he meanwhile i succeeded to compile the thrust prgrams I had to use boost1.461 instead of 1.480
But 1 Example still is crashing:
mma@arnie /downloads/odeintv2/libs/numeric/odeint/examples/thrust> ./phase_oscillator_chain
Segmentation fault
the others seem to run very slow, for example lorentz_parameters: i interrupted after 23minutes:
113
114
115
116
117
^C
real 23m38.396s
user 186m48.651s
sys 0m15.544s
is this normal?
(i really have 2 gtx580 nvidia cards builtin)
lspci  grep i nvidia
01:00.0 PCI bridge: nVidia Corporation Device 05bf (rev a3)
02:00.0 PCI bridge: nVidia Corporation Device 05bf (rev a3)
02:01.0 PCI bridge: nVidia Corporation Device 05bf (rev a3)
02:02.0 PCI bridge: nVidia Corporation Device 05bf (rev a3)
02:03.0 PCI bridge: nVidia Corporation Device 05bf (rev a3)
03:00.0 VGA compatible controller: nVidia Corporation GF110 [GeForce GTX 580] (rev a1)
03:00.1 Audio device: nVidia Corporation GF110 High Definition Audio Controller (rev a1)
05:00.0 VGA compatible controller: nVidia Corporation GF110 [GeForce GTX 580] (rev a1)
05:00.1 Audio device: nVidia Corporation GF110 High Definition Audio Controller (rev a1
)
Greetings Michael





Sorry for the late reply. The execution time is not normal, it should be much faster.
I guess you already compile with full optimization and for the GPU instead for OpenMP. The Makefile distributed with odeint builds for OpenMP.
I dont have a GPU available right now, I have to check later today.





hello
well because i am interested in solving odes on a gpu i was interested to see how this works with "thrust"
as i understand from your article thrust pushed and pulls the code (ordinarily spoken) from and to the gpu
so i wanted to see how this works.
watching the gpu indead show that she is idle most of the time.
but the cpu (i72600 4 cores with 2 threads) show that all 8 virtual cpus are working fullspeed)
so this seems to work. But how can i now test with the gpu?
i compiled with
NVCCFLAGS = O3 $(INCLUDES) arch $(ARCH) compilerbindir=/usr/bin/g++ Xcompiler fopenmp DTHRUST_DEVICE_BACKEND=THRUST_DEVICE_BACKEND_OMP





Ok, this means you are using OpenMP and the long run times seem to be normal. (Thats a nice feature of thrust: you can use it with OpenMP).
To build for the GPU simply remove Xcompiler fopenmp DTHRUST_DEVICE_BACKEND=THRUST_DEVICE_BACKEND_OMP from the NVCCFLAGS. That should be all.





hello
i changed the makefile and tried the phase_oscillator_ensemble and it indeed used the gpu:
.............
.............
.............
Dopri5 : 4.7 0.755809 9.44 643 3546
Dopri5 : 4.8 0.761798 9.55 643 3617
Dopri5 : 4.9 0.767404 9.93 664 3701
Dopri5 : 5 0.772688 10.04 678 3771
Segmentation fault
it needed about 34 minutes. is this normal (besides the segmentation fault)?
regards Michael





How long took the LorenzParameter example?





still running 6 minutes and the screen shows the number 413
(the loop itself in source code is 10.000)
that means about 0.8sec for 1 step or nearly 3 h for the whole loop.
But the gpu is running at 69% speed and its tempareture is 87 degree celsius (=188 dgree fahrenheit)
so the gpu seems to work with the code.
but 0.8s for 1 step is this possible?
(.. in the meanwhile step 950 and 17min44sec running)
greetings michael





So, you have a speedup factor of approx. 10. This looks ok. In the meanwhile we found some bug in thrust_resizer which might be responsible for the crashes of dopri5 in the phase_oscillator_ensemble example. I am not sure if it is fixed right now, but will keep you informed.





nice work. I got the same result with ubuntu linux 10.10.





i just commited a bugfix for the thrust backend in odeint. i hope this helps with your segmentation faults. sorry for this, i forgot to adjust the thrust backend when applying some changes in the memory management.
regards, mario





would you explain the method using cuda and thrust with the simplest mathematical ODEs? because physical relationships in all three examples are confusing for me.





Sorry for the late reply. Solving a very simple ODE does not make sense with CUDA, unless you are solve many of them in parallel. What kind of ODE are you familar with? I could provide an example for dxdt = a x, where several ODEs with different parameters a are solved in parallel. Is this ok for you?





A single unit of my question is as follows:
dydx[1]=dydx[2];
dydx[2]=f(x);
dydx[3]=g(x).
f and g are functions about x. There are many group in my question. but if you provide dydx = ax is ok.
modified 4Dec11 6:22am.






It is enough to provide me a start for understanding all thrust examples. Thank you. there are too many materials to be finished in this month.






download boost_1_47_0.zip and odeintv2.zip
unrar odeintv2.zip to e:\odeintv2，like E:\odeintv2\boost\numeric\odeint.hpp
unrar boost_1_47_0.zip to e:\odeintv2，like E:\odeintv2\boost\boost.png，rename several files with same file name。
place E:\odeintv2 into VS2008 tools – options  c++ Directories  include files
create a empty win32 command project, add a main.cpp, filling with following content:
#include <iostream>
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>
using namespace std;
using namespace boost::numeric::odeint;
const double sigma = 10.0;
const double R = 28.0;
const double b = 8.0 / 3.0;
typedef boost::array< double , 3 > state_type;
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] = b * x[2] + x[0] * x[1];
}
void write_lorenz( const state_type &x , const double t )
{
cout << t << '\t' << x[0] << '\t' << x[1] << '\t' << x[2] << endl;
}
int main(int argc, char **argv)
{
state_type x = { 10.0 , 1.0 , 1.0 }; // initial conditions
integrate( lorenz , x , 0.0 , 25.0 , 0.1 , write_lorenz );
}
compile, ignore warnings and run.







General News Suggestion Question Bug Answer Joke Rant Admin Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

Type  Article 
Licence  CPOL 
First Posted  19 Oct 2011 
Views  40,341 
Bookmarked  30 times 

