ceph/src/boost/libs/numeric/odeint/test_external/gmp/gmp_integrate.cpp

   1 /* Boost check_gmp.cpp test file
   2
   3  Copyright 2010-2012 Mario Mulansky
   4  Copyright 2011-2012 Karsten Ahnert
   5
   6  This file tests the odeint library with the gmp arbitrary precision types
   7
   8  Distributed under the Boost Software License, Version 1.0.
   9  (See accompanying file LICENSE_1_0.txt or
  10  copy at http://www.boost.org/LICENSE_1_0.txt)
  11 */
  12
  13 #define BOOST_TEST_MODULE odeint_gmp
  14
  15 #include <gmpxx.h>
  16
  17 #include <boost/test/unit_test.hpp>
  18 #include <boost/array.hpp>
  19
  20 #include <boost/mpl/vector.hpp>
  21
  22 #include <boost/numeric/odeint.hpp>
  23 //#include <boost/numeric/odeint/algebra/vector_space_algebra.hpp>
  24
  25 using namespace boost::unit_test;
  26 using namespace boost::numeric::odeint;
  27
  28 namespace mpl = boost::mpl;
  29
  30 const int precision = 1024;
  31
  32 typedef mpf_class value_type;
  33 typedef mpf_class state_type;
  34
  35 //provide min, max and pow functions for mpf types - required for controlled steppers
  36 value_type min( const value_type a , const value_type b )
  37 {
  38     if( a<b ) return a;
  39     else return b;
  40 }
  41 value_type max( const value_type a , const value_type b )
  42 {
  43     if( a>b ) return a;
  44     else return b;
  45 }
  46 value_type pow( const value_type a , const value_type b )
  47 {
  48     // do the calculation in double precision...
  49     return value_type( std::pow( a.get_d() , b.get_d() ) );
  50 }
  51
  52
  53 //provide vector_space reduce:
  54
  55 namespace boost { namespace numeric { namespace odeint {
  56
  57 template<>
  58 struct vector_space_reduce< state_type >
  59 {
  60   template< class Op >
  61   state_type operator()( state_type x , Op op , state_type init ) const
  62   {
  63       init = op( init , x );
  64       return init;
  65   }
  66 };
  67
  68 } } }
  69
  70
  71 void constant_system( const state_type &x , state_type &dxdt , value_type t )
  72 {
  73     dxdt = value_type( 1.0 , precision );
  74 }
  75
  76 /* check runge kutta stepers */
  77 typedef mpl::vector<
  78     euler< state_type , value_type , state_type , value_type , vector_space_algebra > ,
  79     modified_midpoint< state_type , value_type , state_type , value_type , vector_space_algebra > ,
  80     runge_kutta4< state_type , value_type , state_type , value_type , vector_space_algebra > ,
  81     runge_kutta4_classic< state_type , value_type , state_type , value_type , vector_space_algebra > ,
  82     runge_kutta_cash_karp54_classic< state_type , value_type , state_type , value_type , vector_space_algebra > ,
  83     runge_kutta_cash_karp54< state_type , value_type , state_type , value_type , vector_space_algebra > ,
  84     runge_kutta_dopri5< state_type , value_type , state_type , value_type , vector_space_algebra > ,
  85     runge_kutta_fehlberg78< state_type , value_type , state_type , value_type , vector_space_algebra >
  86     > stepper_types;
  87
  88
  89 template< class Stepper >
  90 struct perform_integrate_const_test {
  91
  92     void operator()( void )
  93     {
  94         /* We have to specify the desired precision in advance! */
  95         mpf_set_default_prec( precision );
  96
  97         mpf_t eps_ , unity;
  98         mpf_init( eps_ ); mpf_init( unity );
  99         mpf_set_d( unity , 1.0 );
 100         mpf_div_2exp( eps_ , unity , precision-1 ); // 2^(-precision+1) : smallest number that can be represented with used precision
 101         value_type eps( eps_ );
 102
 103         Stepper stepper;
 104         state_type x;
 105         x = 0.0;
 106         value_type t0( 0.0 );
 107         value_type tend( 1.0 );
 108         value_type dt(0.1);
 109
 110         integrate_const( stepper , constant_system , x , t0 , tend , dt );
 111
 112         x = 0.0;
 113         t0 = 0.0;
 114         dt = 0.1;
 115         size_t steps = 10;
 116
 117         integrate_n_steps( stepper , constant_system , x , t0 , dt , steps );
 118
 119         BOOST_CHECK_MESSAGE( abs( x - 10*dt ) < eps , x );
 120     }
 121 };
 122
 123
 124 BOOST_AUTO_TEST_CASE_TEMPLATE( integrate_const_test , Stepper , stepper_types )
 125 {
 126     perform_integrate_const_test< Stepper > tester;
 127     tester();
 128 }
 129
 130
 131 typedef mpl::vector<
 132     controlled_runge_kutta< runge_kutta_cash_karp54_classic< state_type , value_type , state_type , value_type , vector_space_algebra > > ,
 133     controlled_runge_kutta< runge_kutta_dopri5< state_type , value_type , state_type , value_type , vector_space_algebra > > ,
 134     controlled_runge_kutta< runge_kutta_fehlberg78< state_type , value_type , state_type , value_type , vector_space_algebra > > ,
 135     bulirsch_stoer< state_type , value_type , state_type , value_type , vector_space_algebra >
 136     > controlled_stepper_types;
 137
 138
 139 template< class Stepper >
 140 struct perform_integrate_adaptive_test {
 141
 142     void operator()( void )
 143     {
 144         mpf_set_default_prec( precision );
 145
 146         mpf_t eps_ , unity;
 147         mpf_init( eps_ ); mpf_init( unity );
 148         mpf_set_d( unity , 1.0 );
 149         mpf_div_2exp( eps_ , unity , precision-1 ); // 2^(-precision+1) : smallest number that can be represented with used precision
 150         value_type eps( eps_ );
 151
 152         Stepper stepper;
 153         state_type x;
 154         x = 0.0;
 155         value_type t0( 0.0 );
 156         value_type tend( 1.0 );
 157         value_type dt(0.1);
 158
 159         integrate_adaptive( stepper , constant_system , x , t0 , tend , dt );
 160
 161         BOOST_CHECK_MESSAGE( abs( x - tend ) < eps , x - 0.1 );
 162     }
 163 };
 164
 165 BOOST_AUTO_TEST_CASE_TEMPLATE( integrate_adaptive__test , Stepper , controlled_stepper_types )
 166 {
 167     perform_integrate_adaptive_test< Stepper > tester;
 168     tester();
 169 }