ceph/src/boost/libs/numeric/odeint/include/boost/numeric/odeint/external/mtl4/implicit_euler_mtl4.hpp

   1 /*
   2 [begin_description]
   3 Modification of the implicit Euler method, works with the MTL4 matrix library only.
   4 [end_description]
   5
   6 Copyright 2012-2013 Andreas Angelopoulos
   7 Copyright 2012-2013 Karsten Ahnert
   8 Copyright 2012-2013 Mario Mulansky
   9
  10 Distributed under the Boost Software License, Version 1.0.
  11 (See accompanying file LICENSE_1_0.txt or
  12 copy at http://www.boost.org/LICENSE_1_0.txt)
  13 */
  14
  15
  16 #ifndef BOOST_NUMERIC_ODEINT_EXTERNAL_IMPLICIT_EULER_MTL4_HPP_INCLUDED
  17 #define BOOST_NUMERIC_ODEINT_EXTERNAL_IMPLICIT_EULER_MTL4_HPP_INCLUDED
  18
  19
  20 #include <utility>
  21
  22 #include <boost/numeric/odeint/util/bind.hpp>
  23 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
  24 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
  25
  26 #include <boost/numeric/odeint/external/mtl4/mtl4_resize.hpp>
  27
  28 #include <boost/numeric/mtl/mtl.hpp>
  29 #include <boost/numeric/itl/itl.hpp>
  30
  31
  32
  33
  34 namespace boost {
  35 namespace numeric {
  36 namespace odeint {
  37
  38
  39 template< class ValueType , class Resizer = initially_resizer >
  40 class implicit_euler_mtl4
  41 {
  42
  43 public:
  44
  45     typedef ValueType value_type;
  46     typedef value_type time_type;
  47     typedef mtl::dense_vector<value_type> state_type;
  48
  49     typedef state_wrapper< state_type > wrapped_state_type;
  50     typedef state_type deriv_type;
  51     typedef state_wrapper< deriv_type > wrapped_deriv_type;
  52     typedef mtl::compressed2D< value_type > matrix_type;
  53     typedef state_wrapper< matrix_type > wrapped_matrix_type;
  54
  55     typedef Resizer resizer_type;
  56     typedef stepper_tag stepper_category;
  57
  58     typedef implicit_euler_mtl4< ValueType , Resizer > stepper_type;
  59
  60
  61     implicit_euler_mtl4( const value_type epsilon = 1E-6 )
  62         : m_epsilon( epsilon ) , m_resizer() ,
  63           m_dxdt() , m_x() ,
  64           m_identity() , m_jacobi()
  65     { }
  66
  67
  68     template< class System >
  69     void do_step( System system , state_type &x , time_type t , time_type dt )
  70     {
  71         typedef typename odeint::unwrap_reference< System >::type system_type;
  72         typedef typename odeint::unwrap_reference< typename system_type::first_type >::type deriv_func_type;
  73         typedef typename odeint::unwrap_reference< typename system_type::second_type >::type jacobi_func_type;
  74         system_type &sys = system;
  75         deriv_func_type &deriv_func = sys.first;
  76         jacobi_func_type &jacobi_func = sys.second;
  77
  78         m_resizer.adjust_size( x , detail::bind(
  79                                    &stepper_type::template resize_impl< state_type > , detail::ref( *this ) , detail::_1 ) );
  80
  81         m_identity.m_v = 1;
  82
  83         t += dt;
  84         m_x.m_v = x;
  85
  86         deriv_func( x , m_dxdt.m_v , t );
  87         jacobi_func( x , m_jacobi.m_v , t );
  88
  89
  90         m_dxdt.m_v *= -dt;
  91
  92         m_jacobi.m_v *= dt;
  93         m_jacobi.m_v -= m_identity.m_v ;
  94
  95
  96
  97         // using ilu_0 preconditioning -incomplete LU factorisation
  98         // itl::pc::diagonal<matrix_type,double> L(m_jacobi.m_v);
  99         itl::pc::ilu_0<matrix_type> L( m_jacobi.m_v );
 100
 101         solve( m_jacobi.m_v , m_x.m_v , m_dxdt.m_v , L );
 102         x+= m_x.m_v;
 103
 104
 105     }
 106
 107
 108     template< class StateType >
 109     void adjust_size( const StateType &x )
 110     {
 111         resize_impl( x );
 112     }
 113
 114
 115 private:
 116
 117
 118     /*
 119       Applying approximate iterative linear solvers
 120       default solver is Biconjugate gradient stabilized method
 121       itl::bicgstab(A, x, b, L, iter);
 122     */
 123     template < class LinearOperator, class HilbertSpaceX, class HilbertSpaceB, class Preconditioner>
 124     void solve(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
 125                const Preconditioner& L, int max_iteractions =500)
 126     {
 127         // Termination criterion: r < 1e-6 * b or N iterations
 128         itl::basic_iteration< double > iter( b , max_iteractions , 1e-6 );
 129         itl::bicgstab( A , x , b , L , iter );
 130
 131     }
 132
 133
 134     template< class StateIn >
 135     bool resize_impl( const StateIn &x )
 136     {
 137         bool resized = false;
 138         resized |= adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
 139         resized |= adjust_size_by_resizeability( m_x , x , typename is_resizeable<state_type>::type() );
 140         resized |= adjust_size_by_resizeability( m_identity , x , typename is_resizeable<matrix_type>::type() );
 141         resized |= adjust_size_by_resizeability( m_jacobi , x , typename is_resizeable<matrix_type>::type() );
 142         return resized;
 143     }
 144
 145
 146 private:
 147
 148     value_type m_epsilon;
 149     resizer_type m_resizer;
 150     wrapped_deriv_type m_dxdt;
 151     wrapped_state_type m_x;
 152     wrapped_matrix_type m_identity;
 153     wrapped_matrix_type m_jacobi;
 154 };
 155
 156
 157 } // odeint
 158 } // numeric
 159 } // boost
 160
 161
 162 #endif // BOOST_NUMERIC_ODEINT_EXTERNAL_IMPLICIT_EULER_MTL4_HPP_INCLUDED