123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324 |
- /*
- * Copyright 2012 Karsten Ahnert
- * Copyright 2012 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)
- */
- #include <iostream>
- #include <fstream>
- #include <utility>
- #include "time.h"
- #include <boost/numeric/odeint.hpp>
- #include <boost/phoenix/phoenix.hpp>
- #include <boost/numeric/mtl/mtl.hpp>
- #include <boost/numeric/odeint/external/mtl4/implicit_euler_mtl4.hpp>
- using namespace std;
- using namespace boost::numeric::odeint;
- namespace phoenix = boost::phoenix;
- typedef mtl::dense_vector< double > vec_mtl4;
- typedef mtl::compressed2D< double > mat_mtl4;
- typedef boost::numeric::ublas::vector< double > vec_ublas;
- typedef boost::numeric::ublas::matrix< double > mat_ublas;
- // two systems defined 1 & 2 both are mostly sparse with the number of element variable
- struct system1_mtl4
- {
- void operator()( const vec_mtl4 &x , vec_mtl4 &dxdt , double t )
- {
- int size = mtl::size(x);
- dxdt[ 0 ] = -0.06*x[0];
- for (int i =1; i< size ; ++i){
- dxdt[ i ] = 4.2*x[i-1]-2.2*x[i]*x[i];
- }
- }
- };
- struct jacobi1_mtl4
- {
- void operator()( const vec_mtl4 &x , mat_mtl4 &J , const double &t )
- {
- int size = mtl::size(x);
- mtl::matrix::inserter<mat_mtl4> ins(J);
- ins[0][0]=-0.06;
- for (int i =1; i< size ; ++i)
- {
- ins[i][i-1] = + 4.2;
- ins[i][i] = -4.2*x[i];
- }
- }
- };
- struct system1_ublas
- {
- void operator()( const vec_ublas &x , vec_ublas &dxdt , double t )
- {
- int size = x.size();
- dxdt[ 0 ] = -0.06*x[0];
- for (int i =1; i< size ; ++i){
- dxdt[ i ] = 4.2*x[i-1]-2.2*x[i]*x[i];
- }
- }
- };
- struct jacobi1_ublas
- {
- void operator()( const vec_ublas &x , mat_ublas &J , const double &t )
- {
- int size = x.size();
- // mtl::matrix::inserter<mat_mtl4> ins(J);
- J(0,0)=-0.06;
- for (int i =1; i< size ; ++i){
- //ins[i][0]=120.0*x[i];
- J(i,i-1) = + 4.2;
- J(i,i) = -4.2*x[i];
- }
- }
- };
- struct system2_mtl4
- {
- void operator()( const vec_mtl4 &x , vec_mtl4 &dxdt , double t )
- {
- int size = mtl::size(x);
- for (int i =0; i< size/5 ; i+=5){
- dxdt[ i ] = -0.5*x[i];
- dxdt[i+1]= +25*x[i+1]*x[i+2]-740*x[i+3]*x[i+3]+4.2e-2*x[i];
- dxdt[i+2]= +25*x[i]*x[i]-740*x[i+3]*x[i+3];
- dxdt[i+3]= -25*x[i+1]*x[i+2]+740*x[i+3]*x[i+3];
- dxdt[i+4] = 0.250*x[i]*x[i+1]-44.5*x[i+3];
- }
- }
- };
- struct jacobi2_mtl4
- {
- void operator()( const vec_mtl4 &x , mat_mtl4 &J , const double &t )
- {
- int size = mtl::size(x);
- mtl::matrix::inserter<mat_mtl4> ins(J);
- for (int i =0; i< size/5 ; i+=5){
- ins[ i ][i] = -0.5;
- ins[i+1][i+1]=25*x[i+2];
- ins[i+1][i+2] = 25*x[i+1];
- ins[i+1][i+3] = -740*2*x[i+3];
- ins[i+1][i] =+4.2e-2;
- ins[i+2][i]= 50*x[i];
- ins[i+2][i+3]= -740*2*x[i+3];
- ins[i+3][i+1] = -25*x[i+2];
- ins[i+3][i+2] = -25*x[i+1];
- ins[i+3][i+3] = +740*2*x[i+3];
- ins[i+4][i] = 0.25*x[i+1];
- ins[i+4][i+1] =0.25*x[i];
- ins[i+4][i+3]=-44.5;
- }
- }
- };
- struct system2_ublas
- {
- void operator()( const vec_ublas &x , vec_ublas &dxdt , double t )
- {
- int size = x.size();
- for (int i =0; i< size/5 ; i+=5){
- dxdt[ i ] = -4.2e-2*x[i];
- dxdt[i+1]= +25*x[i+1]*x[i+2]-740*x[i+3]*x[i+3]+4.2e-2*x[i];
- dxdt[i+2]= +25*x[i]*x[i]-740*x[i+3]*x[i+3];
- dxdt[i+3]= -25*x[i+1]*x[i+2]+740*x[i+3]*x[i+3];
- dxdt[i+4] = 0.250*x[i]*x[i+1]-44.5*x[i+3];
- }
- }
- };
- struct jacobi2_ublas
- {
- void operator()( const vec_ublas &x , mat_ublas &J , const double &t )
- {
- int size = x.size();
- for (int i =0; i< size/5 ; i+=5){
- J(i ,i) = -4.2e-2;
- J(i+1,i+1)=25*x[i+2];
- J(i+1,i+2) = 25*x[i+1];
- J(i+1,i+3) = -740*2*x[i+3];
- J(i+1,i) =+4.2e-2;
- J(i+2,i)= 50*x[i];
- J(i+2,i+3)= -740*2*x[i+3];
- J(i+3,i+1) = -25*x[i+2];
- J(i+3,i+2) = -25*x[i+1];
- J(i+3,i+3) = +740*2*x[i+3];
- J(i+4,i) = 0.25*x[i+1];
- J(i+4,i+1) =0.25*x[i];
- J(i+4,i+3)=-44.5;
- }
- }
- };
- void testRidiculouslyMassiveArray( int size )
- {
- typedef boost::numeric::odeint::implicit_euler_mtl4 < double > mtl4stepper;
- typedef boost::numeric::odeint::implicit_euler< double > booststepper;
- vec_mtl4 x(size , 0.0);
- x[0]=1;
- double dt = 0.02;
- double endtime = 10.0;
- clock_t tstart_mtl4 = clock();
- size_t num_of_steps_mtl4 = integrate_const(
- mtl4stepper() ,
- make_pair( system1_mtl4() , jacobi1_mtl4() ) ,
- x , 0.0 , endtime , dt );
- clock_t tend_mtl4 = clock() ;
- clog << x[0] << endl;
- clog << num_of_steps_mtl4 << " time elapsed: " << (double)(tend_mtl4-tstart_mtl4 )/CLOCKS_PER_SEC << endl;
- vec_ublas x_ublas(size , 0.0);
- x_ublas[0]=1;
- clock_t tstart_boost = clock();
- size_t num_of_steps_ublas = integrate_const(
- booststepper() ,
- make_pair( system1_ublas() , jacobi1_ublas() ) ,
- x_ublas , 0.0 , endtime , dt );
- clock_t tend_boost = clock() ;
-
- clog << x_ublas[0] << endl;
- clog << num_of_steps_ublas << " time elapsed: " << (double)(tend_boost-tstart_boost)/CLOCKS_PER_SEC<< endl;
- clog << "dt_ublas/dt_mtl4 = " << (double)( tend_boost-tstart_boost )/( tend_mtl4-tstart_mtl4 ) << endl << endl;
- return ;
- }
- void testRidiculouslyMassiveArray2( int size )
- {
- typedef boost::numeric::odeint::implicit_euler_mtl4 < double > mtl4stepper;
- typedef boost::numeric::odeint::implicit_euler< double > booststepper;
- vec_mtl4 x(size , 0.0);
- x[0]=100;
- double dt = 0.01;
- double endtime = 10.0;
- clock_t tstart_mtl4 = clock();
- size_t num_of_steps_mtl4 = integrate_const(
- mtl4stepper() ,
- make_pair( system1_mtl4() , jacobi1_mtl4() ) ,
- x , 0.0 , endtime , dt );
- clock_t tend_mtl4 = clock() ;
-
- clog << x[0] << endl;
- clog << num_of_steps_mtl4 << " time elapsed: " << (double)(tend_mtl4-tstart_mtl4 )/CLOCKS_PER_SEC << endl;
- vec_ublas x_ublas(size , 0.0);
- x_ublas[0]=100;
- clock_t tstart_boost = clock();
- size_t num_of_steps_ublas = integrate_const(
- booststepper() ,
- make_pair( system1_ublas() , jacobi1_ublas() ) ,
- x_ublas , 0.0 , endtime , dt );
- clock_t tend_boost = clock() ;
-
- clog << x_ublas[0] << endl;
- clog << num_of_steps_ublas << " time elapsed: " << (double)(tend_boost-tstart_boost)/CLOCKS_PER_SEC<< endl;
- clog << "dt_ublas/dt_mtl4 = " << (double)( tend_boost-tstart_boost )/( tend_mtl4-tstart_mtl4 ) << endl << endl;
- return ;
- }
-
- int main( int argc , char **argv )
- {
- std::vector< size_t > length;
- length.push_back( 8 );
- length.push_back( 16 );
- length.push_back( 32 );
- length.push_back( 64 );
- length.push_back( 128 );
- length.push_back( 256 );
- for( size_t i=0 ; i<length.size() ; ++i )
- {
- clog << "Testing with size " << length[i] << endl;
- testRidiculouslyMassiveArray( length[i] );
- }
- clog << endl << endl;
- for( size_t i=0 ; i<length.size() ; ++i )
- {
- clog << "Testing with size " << length[i] << endl;
- testRidiculouslyMassiveArray2( length[i] );
- }
- }
|