1 /*
2  [auto_generated]
3  boost/numeric/odeint/util/ublas_wrapper.hpp
4 
5  [begin_description]
6  Resizing for ublas::vector and ublas::matrix
7  [end_description]
8 
9  Copyright 2011-2013 Mario Mulansky
10  Copyright 2011-2013 Karsten Ahnert
11 
12  Distributed under the Boost Software License, Version 1.0.
13  (See accompanying file LICENSE_1_0.txt or
14  copy at http://www.boost.org/LICENSE_1_0.txt)
15  */
16 
17 
18 #ifndef BOOST_NUMERIC_ODEINT_UTIL_UBLAS_WRAPPER_HPP_INCLUDED
19 #define BOOST_NUMERIC_ODEINT_UTIL_UBLAS_WRAPPER_HPP_INCLUDED
20 
21 
22 #include <boost/type_traits/integral_constant.hpp>
23 #include <boost/numeric/ublas/vector.hpp>
24 #include <boost/numeric/ublas/matrix.hpp>
25 #include <boost/numeric/ublas/lu.hpp>
26 #include <boost/numeric/ublas/vector_expression.hpp>
27 #include <boost/numeric/ublas/matrix_expression.hpp>
28 
29 #include <boost/numeric/odeint/algebra/vector_space_algebra.hpp>
30 #include <boost/numeric/odeint/algebra/default_operations.hpp>
31 
32 #include <boost/numeric/odeint/util/is_resizeable.hpp>
33 #include <boost/numeric/odeint/util/state_wrapper.hpp>
34 
35 
36 /* extend ublas by a few operations */
37 
38 /* map norm_inf onto reduce( v , default_operations::maximum ) */
39 namespace boost { namespace numeric { namespace odeint {
40 
41     template< typename T , typename A >
42     struct vector_space_norm_inf< boost::numeric::ublas::vector<T,A> >
43     {
44         typedef T result_type;
45 
operator ()boost::numeric::odeint::vector_space_norm_inf46         result_type operator()( const boost::numeric::ublas::vector<T,A> &x ) const
47         {
48             return boost::numeric::ublas::norm_inf( x );
49         }
50     };
51 
52 
53     template< class T , class L , class A >
54     struct vector_space_norm_inf< boost::numeric::ublas::matrix<T,L,A> >
55     {
56         typedef T result_type;
57 
operator ()boost::numeric::odeint::vector_space_norm_inf58         result_type operator()( const boost::numeric::ublas::matrix<T,L,A> &x ) const
59         {
60             return boost::numeric::ublas::norm_inf( x );
61         }
62     };
63 } } }
64 
65 /* additional operations:
66  * abs( v )
67  * v / w
68  * a + v
69  */
70 namespace boost { namespace numeric { namespace ublas {
71 
72 
73     // elementwise abs - calculates absolute values of the elements
74     template<class T>
75     struct scalar_abs: public scalar_unary_functor<T> {
76         typedef typename scalar_unary_functor<T>::value_type value_type;
77         typedef typename scalar_unary_functor<T>::argument_type argument_type;
78         typedef typename scalar_unary_functor<T>::result_type result_type;
79 
80         static BOOST_UBLAS_INLINE
applyboost::numeric::ublas::scalar_abs81         result_type apply (argument_type t) {
82             using std::abs;
83             return abs (t);
84         }
85     };
86 
87 
88     // (abs v) [i] = abs (v [i])
89     template<class E>
90     BOOST_UBLAS_INLINE
91     typename vector_unary_traits<E, scalar_abs<typename E::value_type> >::result_type
abs(const vector_expression<E> & e)92     abs (const vector_expression<E> &e) {
93         typedef typename vector_unary_traits<E, scalar_abs<typename E::value_type> >::expression_type expression_type;
94         return expression_type (e ());
95     }
96 
97     // (abs m) [i] = abs (m [i])
98     template<class E>
99     BOOST_UBLAS_INLINE
100     typename matrix_unary1_traits<E, scalar_abs<typename E::value_type> >::result_type
abs(const matrix_expression<E> & e)101     abs (const matrix_expression<E> &e) {
102         typedef typename matrix_unary1_traits<E, scalar_abs<typename E::value_type> >::expression_type expression_type;
103         return expression_type (e ());
104     }
105 
106 
107     // elementwise division (v1 / v2) [i] = v1 [i] / v2 [i]
108     template<class E1, class E2>
109     BOOST_UBLAS_INLINE
110     typename vector_binary_traits<E1, E2, scalar_divides<typename E1::value_type,
111                                                          typename E2::value_type> >::result_type
operator /(const vector_expression<E1> & e1,const vector_expression<E2> & e2)112     operator / (const vector_expression<E1> &e1,
113                 const vector_expression<E2> &e2) {
114         typedef typename vector_binary_traits<E1, E2, scalar_divides<typename E1::value_type,
115                                                                      typename E2::value_type> >::expression_type expression_type;
116         return expression_type (e1 (), e2 ());
117     }
118 
119 
120     // elementwise division (m1 / m2) [i] = m1 [i] / m2 [i]
121     template<class E1, class E2>
122     BOOST_UBLAS_INLINE
123     typename matrix_binary_traits<E1, E2, scalar_divides<typename E1::value_type,
124     typename E2::value_type> >::result_type
operator /(const matrix_expression<E1> & e1,const matrix_expression<E2> & e2)125     operator / (const matrix_expression<E1> &e1,
126                 const matrix_expression<E2> &e2) {
127         typedef typename matrix_binary_traits<E1, E2, scalar_divides<typename E1::value_type,
128                                                                      typename E2::value_type> >::expression_type expression_type;
129         return expression_type (e1 (), e2 ());
130     }
131 
132     // addition with scalar
133     // (t + v) [i] = t + v [i]
134     template<class T1, class E2>
135     BOOST_UBLAS_INLINE
136     typename enable_if< is_convertible<T1, typename E2::value_type >,
137     typename vector_binary_scalar1_traits<const T1, E2, scalar_plus<T1, typename E2::value_type> >::result_type
138     >::type
operator +(const T1 & e1,const vector_expression<E2> & e2)139     operator + (const T1 &e1,
140                 const vector_expression<E2> &e2) {
141         typedef typename vector_binary_scalar1_traits<const T1, E2, scalar_plus<T1, typename E2::value_type> >::expression_type expression_type;
142         return expression_type (e1, e2 ());
143     }
144 
145     // addition with scalar
146     // (t + m) [i] = t + m [i]
147     template<class T1, class E2>
148     BOOST_UBLAS_INLINE
149     typename enable_if< is_convertible<T1, typename E2::value_type >,
150     typename matrix_binary_scalar1_traits<const T1, E2, scalar_plus<T1, typename E2::value_type> >::result_type
151     >::type
operator +(const T1 & e1,const matrix_expression<E2> & e2)152     operator + (const T1 &e1,
153                 const matrix_expression<E2> &e2) {
154         typedef typename matrix_binary_scalar1_traits<const T1, E2, scalar_plus<T1, typename E2::value_type> >::expression_type expression_type;
155         return expression_type (e1, e2 ());
156     }
157 
158 } } }
159 
160 
161 
162 
163 /* add resize functionality */
164 namespace boost {
165 namespace numeric {
166 namespace odeint {
167 
168 /*
169  * resizeable specialization for boost::numeric::ublas::vector
170  */
171 template< class T , class A >
172 struct is_resizeable< boost::numeric::ublas::vector< T , A > >
173 {
174     typedef boost::true_type type;
175     const static bool value = type::value;
176 };
177 
178 
179 /*
180  * resizeable specialization for boost::numeric::ublas::matrix
181  */
182 template< class T , class L , class A >
183 struct is_resizeable< boost::numeric::ublas::matrix< T , L , A > >
184 {
185     typedef boost::true_type type;
186     const static bool value = type::value;
187 };
188 
189 
190 /*
191  * resizeable specialization for boost::numeric::ublas::permutation_matrix
192  */
193 template< class T , class A >
194 struct is_resizeable< boost::numeric::ublas::permutation_matrix< T , A > >
195 {
196     typedef boost::true_type type;
197     const static bool value = type::value;
198 };
199 
200 
201 // specialization for ublas::matrix
202 // same size and resize specialization for matrix-matrix resizing
203 template< class T , class L , class A , class T2 , class L2 , class A2 >
204 struct same_size_impl< boost::numeric::ublas::matrix< T , L , A > , boost::numeric::ublas::matrix< T2 , L2 , A2 > >
205 {
same_sizeboost::numeric::odeint::same_size_impl206     static bool same_size( const boost::numeric::ublas::matrix< T , L , A > &m1 ,
207                            const boost::numeric::ublas::matrix< T2 , L2 , A2 > &m2 )
208     {
209         return ( ( m1.size1() == m2.size1() ) && ( m1.size2() == m2.size2() ) );
210     }
211 };
212 
213 template< class T , class L , class A , class T2 , class L2 , class A2 >
214 struct resize_impl< boost::numeric::ublas::matrix< T , L , A > , boost::numeric::ublas::matrix< T2 , L2 , A2 > >
215 {
resizeboost::numeric::odeint::resize_impl216     static void resize( boost::numeric::ublas::matrix< T , L , A > &m1 ,
217                         const boost::numeric::ublas::matrix< T2 , L2 , A2 > &m2 )
218     {
219         m1.resize( m2.size1() , m2.size2() );
220     }
221 };
222 
223 
224 
225 // same size and resize specialization for matrix-vector resizing
226 template< class T , class L , class A , class T_V , class A_V >
227 struct same_size_impl< boost::numeric::ublas::matrix< T , L , A > , boost::numeric::ublas::vector< T_V , A_V > >
228 {
same_sizeboost::numeric::odeint::same_size_impl229     static bool same_size( const boost::numeric::ublas::matrix< T , L , A > &m ,
230                            const boost::numeric::ublas::vector< T_V , A_V > &v )
231     {
232         return ( ( m.size1() == v.size() ) && ( m.size2() == v.size() ) );
233     }
234 };
235 
236 template< class T , class L , class A , class T_V , class A_V >
237 struct resize_impl< boost::numeric::ublas::matrix< T , L , A > , boost::numeric::ublas::vector< T_V , A_V > >
238 {
resizeboost::numeric::odeint::resize_impl239     static void resize( boost::numeric::ublas::matrix< T , L , A > &m ,
240                         const boost::numeric::ublas::vector< T_V , A_V > &v )
241     {
242         m.resize( v.size() , v.size() );
243     }
244 };
245 
246 
247 
248 // specialization for ublas::permutation_matrix
249 // same size and resize specialization for matrix-vector resizing
250 template< class T , class A , class T_V , class A_V >
251 struct same_size_impl< boost::numeric::ublas::permutation_matrix< T , A > ,
252                        boost::numeric::ublas::vector< T_V , A_V > >
253 {
same_sizeboost::numeric::odeint::same_size_impl254     static bool same_size( const boost::numeric::ublas::permutation_matrix< T , A > &m ,
255                            const boost::numeric::ublas::vector< T_V , A_V > &v )
256     {
257         return ( m.size() == v.size() ); // && ( m.size2() == v.size() ) );
258     }
259 };
260 
261 template< class T , class A , class T_V , class A_V >
262 struct resize_impl< boost::numeric::ublas::vector< T_V , A_V > ,
263                     boost::numeric::ublas::permutation_matrix< T , A > >
264 {
resizeboost::numeric::odeint::resize_impl265     static void resize( const boost::numeric::ublas::vector< T_V , A_V > &v,
266                         boost::numeric::ublas::permutation_matrix< T , A > &m )
267     {
268         m.resize( v.size() , v.size() );
269     }
270 };
271 
272 
273 
274 
275 
276 
277 
278 template< class T , class A >
279 struct state_wrapper< boost::numeric::ublas::permutation_matrix< T , A > > // with resizing
280 {
281     typedef boost::numeric::ublas::permutation_matrix< T , A > state_type;
282     typedef state_wrapper< state_type > state_wrapper_type;
283 
284     state_type m_v;
285 
state_wrapperboost::numeric::odeint::state_wrapper286     state_wrapper() : m_v( 1 ) // permutation matrix constructor requires a size, choose 1 as default
287     { }
288 
289 };
290 
291 
292 
293 
294 } } }
295 
296 
297 #endif // BOOST_NUMERIC_ODEINT_UTIL_UBLAS_WRAPPER_HPP_INCLUDED
298