odeint/util/ublas_wrapper.hpp

0001 /*
0002  [auto_generated]
0003  boost/numeric/odeint/util/ublas_wrapper.hpp
0004
0005  [begin_description]
0006  Resizing for ublas::vector and ublas::matrix
0007  [end_description]
0008
0009  Copyright 2011-2013 Mario Mulansky
0010  Copyright 2011-2013 Karsten Ahnert
0011
0012  Distributed under the Boost Software License, Version 1.0.
0013  (See accompanying file LICENSE_1_0.txt or
0014  copy at http://www.boost.org/LICENSE_1_0.txt)
0015  */
0016
0017
0018 #ifndef BOOST_NUMERIC_ODEINT_UTIL_UBLAS_WRAPPER_HPP_INCLUDED
0019 #define BOOST_NUMERIC_ODEINT_UTIL_UBLAS_WRAPPER_HPP_INCLUDED
0020
0021
0022 #include <boost/type_traits/integral_constant.hpp>
0023 #include <boost/numeric/ublas/vector.hpp>
0024 #include <boost/numeric/ublas/matrix.hpp>
0025 #include <boost/numeric/ublas/lu.hpp>
0026 #include <boost/numeric/ublas/vector_expression.hpp>
0027 #include <boost/numeric/ublas/matrix_expression.hpp>
0028
0029 #include <boost/numeric/odeint/algebra/vector_space_algebra.hpp>
0030 #include <boost/numeric/odeint/algebra/default_operations.hpp>
0031
0032 #include <boost/numeric/odeint/util/is_resizeable.hpp>
0033 #include <boost/numeric/odeint/util/state_wrapper.hpp>
0034
0035
0036 /* extend ublas by a few operations */
0037
0038 /* map norm_inf onto reduce( v , default_operations::maximum ) */
0039 namespace boost { namespace numeric { namespace odeint {
0040
0041     template< typename T , typename A >
0042     struct vector_space_norm_inf< boost::numeric::ublas::vector<T,A> >
0043     {
0044         typedef T result_type;
0045
0046         result_type operator()( const boost::numeric::ublas::vector<T,A> &x ) const
0047         {
0048             return boost::numeric::ublas::norm_inf( x );
0049         }
0050     };
0051
0052
0053     template< class T , class L , class A >
0054     struct vector_space_norm_inf< boost::numeric::ublas::matrix<T,L,A> >
0055     {
0056         typedef T result_type;
0057
0058         result_type operator()( const boost::numeric::ublas::matrix<T,L,A> &x ) const
0059         {
0060             return boost::numeric::ublas::norm_inf( x );
0061         }
0062     };
0063 } } }
0064
0065 /* additional operations:
0066  * abs( v )
0067  * v / w
0068  * a + v
0069  */
0070 namespace boost { namespace numeric { namespace ublas {
0071
0072
0073     // elementwise abs - calculates absolute values of the elements
0074     template<class T>
0075     struct scalar_abs: public scalar_unary_functor<T> {
0076         typedef typename scalar_unary_functor<T>::value_type value_type;
0077         typedef typename scalar_unary_functor<T>::argument_type argument_type;
0078         typedef typename scalar_unary_functor<T>::result_type result_type;
0079
0080         static BOOST_UBLAS_INLINE
0081         result_type apply (argument_type t) {
0082             using std::abs;
0083             return abs (t);
0084         }
0085     };
0086
0087
0088     // (abs v) [i] = abs (v [i])
0089     template<class E>
0090     BOOST_UBLAS_INLINE
0091     typename vector_unary_traits<E, scalar_abs<typename E::value_type> >::result_type
0092     abs (const vector_expression<E> &e) {
0093         typedef typename vector_unary_traits<E, scalar_abs<typename E::value_type> >::expression_type expression_type;
0094         return expression_type (e ());
0095     }
0096
0097     // (abs m) [i] = abs (m [i])
0098     template<class E>
0099     BOOST_UBLAS_INLINE
0100     typename matrix_unary1_traits<E, scalar_abs<typename E::value_type> >::result_type
0101     abs (const matrix_expression<E> &e) {
0102         typedef typename matrix_unary1_traits<E, scalar_abs<typename E::value_type> >::expression_type expression_type;
0103         return expression_type (e ());
0104     }
0105
0106
0107     // elementwise division (v1 / v2) [i] = v1 [i] / v2 [i]
0108     template<class E1, class E2>
0109     BOOST_UBLAS_INLINE
0110     typename vector_binary_traits<E1, E2, scalar_divides<typename E1::value_type,
0111                                                          typename E2::value_type> >::result_type
0112     operator / (const vector_expression<E1> &e1,
0113                 const vector_expression<E2> &e2) {
0114         typedef typename vector_binary_traits<E1, E2, scalar_divides<typename E1::value_type,
0115                                                                      typename E2::value_type> >::expression_type expression_type;
0116         return expression_type (e1 (), e2 ());
0117     }
0118
0119
0120     // elementwise division (m1 / m2) [i] = m1 [i] / m2 [i]
0121     template<class E1, class E2>
0122     BOOST_UBLAS_INLINE
0123     typename matrix_binary_traits<E1, E2, scalar_divides<typename E1::value_type,
0124     typename E2::value_type> >::result_type
0125     operator / (const matrix_expression<E1> &e1,
0126                 const matrix_expression<E2> &e2) {
0127         typedef typename matrix_binary_traits<E1, E2, scalar_divides<typename E1::value_type,
0128                                                                      typename E2::value_type> >::expression_type expression_type;
0129         return expression_type (e1 (), e2 ());
0130     }
0131
0132     // addition with scalar
0133     // (t + v) [i] = t + v [i]
0134     template<class T1, class E2>
0135     BOOST_UBLAS_INLINE
0136     typename enable_if< is_convertible<T1, typename E2::value_type >,
0137     typename vector_binary_scalar1_traits<const T1, E2, scalar_plus<T1, typename E2::value_type> >::result_type
0138     >::type
0139     operator + (const T1 &e1,
0140                 const vector_expression<E2> &e2) {
0141         typedef typename vector_binary_scalar1_traits<const T1, E2, scalar_plus<T1, typename E2::value_type> >::expression_type expression_type;
0142         return expression_type (e1, e2 ());
0143     }
0144
0145     // addition with scalar
0146     // (t + m) [i] = t + m [i]
0147     template<class T1, class E2>
0148     BOOST_UBLAS_INLINE
0149     typename enable_if< is_convertible<T1, typename E2::value_type >,
0150     typename matrix_binary_scalar1_traits<const T1, E2, scalar_plus<T1, typename E2::value_type> >::result_type
0151     >::type
0152     operator + (const T1 &e1,
0153                 const matrix_expression<E2> &e2) {
0154         typedef typename matrix_binary_scalar1_traits<const T1, E2, scalar_plus<T1, typename E2::value_type> >::expression_type expression_type;
0155         return expression_type (e1, e2 ());
0156     }
0157
0158 } } }
0159
0160
0161
0162
0163 /* add resize functionality */
0164 namespace boost {
0165 namespace numeric {
0166 namespace odeint {
0167
0168 /*
0169  * resizeable specialization for boost::numeric::ublas::vector
0170  */
0171 template< class T , class A >
0172 struct is_resizeable< boost::numeric::ublas::vector< T , A > >
0173 {
0174     typedef boost::true_type type;
0175     const static bool value = type::value;
0176 };
0177
0178
0179 /*
0180  * resizeable specialization for boost::numeric::ublas::matrix
0181  */
0182 template< class T , class L , class A >
0183 struct is_resizeable< boost::numeric::ublas::matrix< T , L , A > >
0184 {
0185     typedef boost::true_type type;
0186     const static bool value = type::value;
0187 };
0188
0189
0190 /*
0191  * resizeable specialization for boost::numeric::ublas::permutation_matrix
0192  */
0193 template< class T , class A >
0194 struct is_resizeable< boost::numeric::ublas::permutation_matrix< T , A > >
0195 {
0196     typedef boost::true_type type;
0197     const static bool value = type::value;
0198 };
0199
0200
0201 // specialization for ublas::matrix
0202 // same size and resize specialization for matrix-matrix resizing
0203 template< class T , class L , class A , class T2 , class L2 , class A2 >
0204 struct same_size_impl< boost::numeric::ublas::matrix< T , L , A > , boost::numeric::ublas::matrix< T2 , L2 , A2 > >
0205 {
0206     static bool same_size( const boost::numeric::ublas::matrix< T , L , A > &m1 ,
0207                            const boost::numeric::ublas::matrix< T2 , L2 , A2 > &m2 )
0208     {
0209         return ( ( m1.size1() == m2.size1() ) && ( m1.size2() == m2.size2() ) );
0210     }
0211 };
0212
0213 template< class T , class L , class A , class T2 , class L2 , class A2 >
0214 struct resize_impl< boost::numeric::ublas::matrix< T , L , A > , boost::numeric::ublas::matrix< T2 , L2 , A2 > >
0215 {
0216     static void resize( boost::numeric::ublas::matrix< T , L , A > &m1 ,
0217                         const boost::numeric::ublas::matrix< T2 , L2 , A2 > &m2 )
0218     {
0219         m1.resize( m2.size1() , m2.size2() );
0220     }
0221 };
0222
0223
0224
0225 // same size and resize specialization for matrix-vector resizing
0226 template< class T , class L , class A , class T_V , class A_V >
0227 struct same_size_impl< boost::numeric::ublas::matrix< T , L , A > , boost::numeric::ublas::vector< T_V , A_V > >
0228 {
0229     static bool same_size( const boost::numeric::ublas::matrix< T , L , A > &m ,
0230                            const boost::numeric::ublas::vector< T_V , A_V > &v )
0231     {
0232         return ( ( m.size1() == v.size() ) && ( m.size2() == v.size() ) );
0233     }
0234 };
0235
0236 template< class T , class L , class A , class T_V , class A_V >
0237 struct resize_impl< boost::numeric::ublas::matrix< T , L , A > , boost::numeric::ublas::vector< T_V , A_V > >
0238 {
0239     static void resize( boost::numeric::ublas::matrix< T , L , A > &m ,
0240                         const boost::numeric::ublas::vector< T_V , A_V > &v )
0241     {
0242         m.resize( v.size() , v.size() );
0243     }
0244 };
0245
0246
0247
0248 // specialization for ublas::permutation_matrix
0249 // same size and resize specialization for matrix-vector resizing
0250 template< class T , class A , class T_V , class A_V >
0251 struct same_size_impl< boost::numeric::ublas::permutation_matrix< T , A > ,
0252                        boost::numeric::ublas::vector< T_V , A_V > >
0253 {
0254     static bool same_size( const boost::numeric::ublas::permutation_matrix< T , A > &m ,
0255                            const boost::numeric::ublas::vector< T_V , A_V > &v )
0256     {
0257         return ( m.size() == v.size() ); // && ( m.size2() == v.size() ) );
0258     }
0259 };
0260
0261 template< class T , class A , class T_V , class A_V >
0262 struct resize_impl< boost::numeric::ublas::vector< T_V , A_V > ,
0263                     boost::numeric::ublas::permutation_matrix< T , A > >
0264 {
0265     static void resize( const boost::numeric::ublas::vector< T_V , A_V > &v,
0266                         boost::numeric::ublas::permutation_matrix< T , A > &m )
0267     {
0268         m.resize( v.size() , v.size() );
0269     }
0270 };
0271
0272
0273
0274
0275
0276
0277
0278 template< class T , class A >
0279 struct state_wrapper< boost::numeric::ublas::permutation_matrix< T , A > > // with resizing
0280 {
0281     typedef boost::numeric::ublas::permutation_matrix< T , A > state_type;
0282     typedef state_wrapper< state_type > state_wrapper_type;
0283
0284     state_type m_v;
0285
0286     state_wrapper() : m_v( 1 ) // permutation matrix constructor requires a size, choose 1 as default
0287     { }
0288
0289 };
0290
0291
0292
0293
0294 } } }
0295
0296
0297 #endif // BOOST_NUMERIC_ODEINT_UTIL_UBLAS_WRAPPER_HPP_INCLUDED