diff options
Diffstat (limited to 'src/boost/libs/numeric/odeint/examples/mtl')
3 files changed, 482 insertions, 0 deletions
diff --git a/src/boost/libs/numeric/odeint/examples/mtl/Jamfile.v2 b/src/boost/libs/numeric/odeint/examples/mtl/Jamfile.v2 new file mode 100644 index 00000000..9f4ffb0f --- /dev/null +++ b/src/boost/libs/numeric/odeint/examples/mtl/Jamfile.v2 @@ -0,0 +1,17 @@ +# Copyright 2011-2013 Mario Mulansky +# Copyright 2012 Karsten Ahnert +# 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) + +# set your MTL4 directory here +MTL4_INCLUDE = /home/mario/MTL4 ; + +project + : requirements + <include>$(MTL4_INCLUDE) + <define>BOOST_ALL_NO_LIB=1 + ; + +exe gauss_packet : gauss_packet.cpp ; +exe implicit_euler_mtl : implicit_euler_mtl.cpp ;
\ No newline at end of file diff --git a/src/boost/libs/numeric/odeint/examples/mtl/gauss_packet.cpp b/src/boost/libs/numeric/odeint/examples/mtl/gauss_packet.cpp new file mode 100644 index 00000000..9579e5d4 --- /dev/null +++ b/src/boost/libs/numeric/odeint/examples/mtl/gauss_packet.cpp @@ -0,0 +1,141 @@ +/* + * gauss_packet.cpp + * + * Schroedinger equation with potential barrier and periodic boundary conditions + * Initial Gauss packet moving to the right + * + * pipe output into gnuplot to see animation + * + * Implementation of Hamilton operator via MTL library + * + * Copyright 2011-2013 Mario Mulansky + * Copyright 2011-2012 Karsten Ahnert + * + * 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 <complex> + +#include <boost/numeric/odeint.hpp> +#include <boost/numeric/odeint/external/mtl4/mtl4.hpp> + +#include <boost/numeric/mtl/mtl.hpp> + + +using namespace std; +using namespace boost::numeric::odeint; + +typedef mtl::dense_vector< complex< double > > state_type; + +struct hamiltonian { + + typedef mtl::compressed2D< complex< double > > matrix_type; + matrix_type m_H; + + hamiltonian( const int N ) : m_H( N , N ) + { + // constructor with zero potential + m_H = 0.0; + initialize_kinetic_term(); + } + + //template< mtl::compressed2D< double > > + hamiltonian( mtl::compressed2D< double > &V ) : m_H( num_rows( V ) , num_cols( V ) ) + { + // use potential V in hamiltonian + m_H = complex<double>( 0.0 , -1.0 ) * V; + initialize_kinetic_term(); + } + + void initialize_kinetic_term( ) + { + const int N = num_rows( m_H ); + mtl::matrix::inserter< matrix_type , mtl::update_plus< complex<double> > > ins( m_H ); + const double z = 1.0; + // fill diagonal and upper and lower diagonal + for( int i = 0 ; i<N ; ++i ) + { + ins[ i ][ (i+1) % N ] << complex< double >( 0.0 , -z ); + ins[ i ][ i ] << complex< double >( 0.0 , z ); + ins[ (i+1) % N ][ i ] << complex< double >( 0.0 , -z ); + } + } + + void operator()( const state_type &psi , state_type &dpsidt , const double t ) + { + dpsidt = m_H * psi; + } + +}; + +struct write_for_gnuplot +{ + size_t m_every , m_count; + + write_for_gnuplot( size_t every = 10 ) + : m_every( every ) , m_count( 0 ) { } + + void operator()( const state_type &x , double t ) + { + if( ( m_count % m_every ) == 0 ) + { + //clog << t << endl; + cout << "p [0:" << mtl::size(x) << "][0:0.02] '-'" << endl; + for( size_t i=0 ; i<mtl::size(x) ; ++i ) + { + cout << i << "\t" << norm(x[i]) << "\n"; + } + cout << "e" << endl; + } + + ++m_count; + } +}; + +static const int N = 1024; +static const int N0 = 256; +static const double sigma0 = 20; +static const double k0 = -1.0; + +int main( int argc , char** argv ) +{ + state_type x( N , 0.0 ); + + // initialize gauss packet with nonzero velocity + for( int i=0 ; i<N ; ++i ) + { + x[i] = exp( -(i-N0)*(i-N0) / ( 4.0*sigma0*sigma0 ) ) * exp( complex< double >( 0.0 , k0*i ) ); + //x[i] += 2.0*exp( -(i+N0-N)*(i+N0-N) / ( 4.0*sigma0*sigma0 ) ) * exp( complex< double >( 0.0 , -k0*i ) ); + } + x /= mtl::two_norm( x ); + + typedef runge_kutta4< state_type > stepper; + + // create potential barrier + mtl::compressed2D< double > V( N , N ); + V = 0.0; + { + mtl::matrix::inserter< mtl::compressed2D< double > > ins( V ); + for( int i=0 ; i<N ; ++i ) + { + //ins[i][i] << 1E-4*(i-N/2)*(i-N/2); + + if( i < N/2 ) + ins[ i ][ i ] << 0.0 ; + else + ins[ i ][ i ] << 1.0 ; + + } + } + + // perform integration, output can be piped to gnuplot + integrate_const( stepper() , hamiltonian( V ) , x , 0.0 , 1000.0 , 0.1 , write_for_gnuplot( 10 ) ); + + clog << "Norm: " << mtl::two_norm( x ) << endl; + + return 0; +} diff --git a/src/boost/libs/numeric/odeint/examples/mtl/implicit_euler_mtl.cpp b/src/boost/libs/numeric/odeint/examples/mtl/implicit_euler_mtl.cpp new file mode 100644 index 00000000..6a1b8e73 --- /dev/null +++ b/src/boost/libs/numeric/odeint/examples/mtl/implicit_euler_mtl.cpp @@ -0,0 +1,324 @@ +/* + * 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] ); + } +} |