1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/os/ossrv/ossrv_pub/boost_apis/boost/numeric/ublas/operation.hpp Fri Jun 15 03:10:57 2012 +0200
1.3 @@ -0,0 +1,855 @@
1.4 +//
1.5 +// Copyright (c) 2000-2002
1.6 +// Joerg Walter, Mathias Koch
1.7 +//
1.8 +// Permission to use, copy, modify, distribute and sell this software
1.9 +// and its documentation for any purpose is hereby granted without fee,
1.10 +// provided that the above copyright notice appear in all copies and
1.11 +// that both that copyright notice and this permission notice appear
1.12 +// in supporting documentation. The authors make no representations
1.13 +// about the suitability of this software for any purpose.
1.14 +// It is provided "as is" without express or implied warranty.
1.15 +//
1.16 +// The authors gratefully acknowledge the support of
1.17 +// GeNeSys mbH & Co. KG in producing this work.
1.18 +//
1.19 +
1.20 +#ifndef _BOOST_UBLAS_OPERATION_
1.21 +#define _BOOST_UBLAS_OPERATION_
1.22 +
1.23 +#include <boost/numeric/ublas/matrix_proxy.hpp>
1.24 +
1.25 +/** \file operation.hpp
1.26 + * \brief This file contains some specialized products.
1.27 + */
1.28 +
1.29 +// axpy-based products
1.30 +// Alexei Novakov had a lot of ideas to improve these. Thanks.
1.31 +// Hendrik Kueck proposed some new kernel. Thanks again.
1.32 +
1.33 +namespace boost { namespace numeric { namespace ublas {
1.34 +
1.35 + template<class V, class T1, class IA1, class TA1, class E2>
1.36 + BOOST_UBLAS_INLINE
1.37 + V &
1.38 + axpy_prod (const compressed_matrix<T1, row_major, 0, IA1, TA1> &e1,
1.39 + const vector_expression<E2> &e2,
1.40 + V &v, row_major_tag) {
1.41 + typedef typename V::size_type size_type;
1.42 + typedef typename V::value_type value_type;
1.43 +
1.44 + for (size_type i = 0; i < e1.filled1 () -1; ++ i) {
1.45 + size_type begin = e1.index1_data () [i];
1.46 + size_type end = e1.index1_data () [i + 1];
1.47 + value_type t (v (i));
1.48 + for (size_type j = begin; j < end; ++ j)
1.49 + t += e1.value_data () [j] * e2 () (e1.index2_data () [j]);
1.50 + v (i) = t;
1.51 + }
1.52 + return v;
1.53 + }
1.54 +
1.55 + template<class V, class T1, class IA1, class TA1, class E2>
1.56 + BOOST_UBLAS_INLINE
1.57 + V &
1.58 + axpy_prod (const compressed_matrix<T1, column_major, 0, IA1, TA1> &e1,
1.59 + const vector_expression<E2> &e2,
1.60 + V &v, column_major_tag) {
1.61 + typedef typename V::size_type size_type;
1.62 +
1.63 + for (size_type j = 0; j < e1.filled1 () -1; ++ j) {
1.64 + size_type begin = e1.index1_data () [j];
1.65 + size_type end = e1.index1_data () [j + 1];
1.66 + for (size_type i = begin; i < end; ++ i)
1.67 + v (e1.index2_data () [i]) += e1.value_data () [i] * e2 () (j);
1.68 + }
1.69 + return v;
1.70 + }
1.71 +
1.72 + // Dispatcher
1.73 + template<class V, class T1, class L1, class IA1, class TA1, class E2>
1.74 + BOOST_UBLAS_INLINE
1.75 + V &
1.76 + axpy_prod (const compressed_matrix<T1, L1, 0, IA1, TA1> &e1,
1.77 + const vector_expression<E2> &e2,
1.78 + V &v, bool init = true) {
1.79 + typedef typename V::value_type value_type;
1.80 + typedef typename L1::orientation_category orientation_category;
1.81 +
1.82 + if (init)
1.83 + v.assign (zero_vector<value_type> (e1.size1 ()));
1.84 +#if BOOST_UBLAS_TYPE_CHECK
1.85 + vector<value_type> cv (v);
1.86 + typedef typename type_traits<value_type>::real_type real_type;
1.87 + real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
1.88 + indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));
1.89 +#endif
1.90 + axpy_prod (e1, e2, v, orientation_category ());
1.91 +#if BOOST_UBLAS_TYPE_CHECK
1.92 + BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
1.93 +#endif
1.94 + return v;
1.95 + }
1.96 + template<class V, class T1, class L1, class IA1, class TA1, class E2>
1.97 + BOOST_UBLAS_INLINE
1.98 + V
1.99 + axpy_prod (const compressed_matrix<T1, L1, 0, IA1, TA1> &e1,
1.100 + const vector_expression<E2> &e2) {
1.101 + typedef V vector_type;
1.102 +
1.103 + vector_type v (e1.size1 ());
1.104 + return axpy_prod (e1, e2, v, true);
1.105 + }
1.106 +
1.107 + template<class V, class T1, class L1, class IA1, class TA1, class E2>
1.108 + BOOST_UBLAS_INLINE
1.109 + V &
1.110 + axpy_prod (const coordinate_matrix<T1, L1, 0, IA1, TA1> &e1,
1.111 + const vector_expression<E2> &e2,
1.112 + V &v, bool init = true) {
1.113 + typedef typename V::size_type size_type;
1.114 + typedef typename V::value_type value_type;
1.115 + typedef L1 layout_type;
1.116 +
1.117 + size_type size1 = e1.size1();
1.118 + size_type size2 = e1.size2();
1.119 +
1.120 + if (init) {
1.121 + noalias(v) = zero_vector<value_type>(size1);
1.122 + }
1.123 +
1.124 + for (size_type i = 0; i < e1.nnz(); ++i) {
1.125 + size_type row_index = layout_type::element1( e1.index1_data () [i], size1, e1.index2_data () [i], size2 );
1.126 + size_type col_index = layout_type::element2( e1.index1_data () [i], size1, e1.index2_data () [i], size2 );
1.127 + v( row_index ) += e1.value_data () [i] * e2 () (col_index);
1.128 + }
1.129 + return v;
1.130 + }
1.131 +
1.132 + template<class V, class E1, class E2>
1.133 + BOOST_UBLAS_INLINE
1.134 + V &
1.135 + axpy_prod (const matrix_expression<E1> &e1,
1.136 + const vector_expression<E2> &e2,
1.137 + V &v, packed_random_access_iterator_tag, row_major_tag) {
1.138 + typedef const E1 expression1_type;
1.139 + typedef const E2 expression2_type;
1.140 + typedef typename V::size_type size_type;
1.141 +
1.142 + typename expression1_type::const_iterator1 it1 (e1 ().begin1 ());
1.143 + typename expression1_type::const_iterator1 it1_end (e1 ().end1 ());
1.144 + while (it1 != it1_end) {
1.145 + size_type index1 (it1.index1 ());
1.146 +#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1.147 + typename expression1_type::const_iterator2 it2 (it1.begin ());
1.148 + typename expression1_type::const_iterator2 it2_end (it1.end ());
1.149 +#else
1.150 + typename expression1_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));
1.151 + typename expression1_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));
1.152 +#endif
1.153 + while (it2 != it2_end) {
1.154 + v (index1) += *it2 * e2 () (it2.index2 ());
1.155 + ++ it2;
1.156 + }
1.157 + ++ it1;
1.158 + }
1.159 + return v;
1.160 + }
1.161 +
1.162 + template<class V, class E1, class E2>
1.163 + BOOST_UBLAS_INLINE
1.164 + V &
1.165 + axpy_prod (const matrix_expression<E1> &e1,
1.166 + const vector_expression<E2> &e2,
1.167 + V &v, packed_random_access_iterator_tag, column_major_tag) {
1.168 + typedef const E1 expression1_type;
1.169 + typedef const E2 expression2_type;
1.170 + typedef typename V::size_type size_type;
1.171 +
1.172 + typename expression1_type::const_iterator2 it2 (e1 ().begin2 ());
1.173 + typename expression1_type::const_iterator2 it2_end (e1 ().end2 ());
1.174 + while (it2 != it2_end) {
1.175 + size_type index2 (it2.index2 ());
1.176 +#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1.177 + typename expression1_type::const_iterator1 it1 (it2.begin ());
1.178 + typename expression1_type::const_iterator1 it1_end (it2.end ());
1.179 +#else
1.180 + typename expression1_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));
1.181 + typename expression1_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));
1.182 +#endif
1.183 + while (it1 != it1_end) {
1.184 + v (it1.index1 ()) += *it1 * e2 () (index2);
1.185 + ++ it1;
1.186 + }
1.187 + ++ it2;
1.188 + }
1.189 + return v;
1.190 + }
1.191 +
1.192 + template<class V, class E1, class E2>
1.193 + BOOST_UBLAS_INLINE
1.194 + V &
1.195 + axpy_prod (const matrix_expression<E1> &e1,
1.196 + const vector_expression<E2> &e2,
1.197 + V &v, sparse_bidirectional_iterator_tag) {
1.198 + typedef const E1 expression1_type;
1.199 + typedef const E2 expression2_type;
1.200 + typedef typename V::size_type size_type;
1.201 +
1.202 + typename expression2_type::const_iterator it (e2 ().begin ());
1.203 + typename expression2_type::const_iterator it_end (e2 ().end ());
1.204 + while (it != it_end) {
1.205 + v.plus_assign (column (e1 (), it.index ()) * *it);
1.206 + ++ it;
1.207 + }
1.208 + return v;
1.209 + }
1.210 +
1.211 + // Dispatcher
1.212 + template<class V, class E1, class E2>
1.213 + BOOST_UBLAS_INLINE
1.214 + V &
1.215 + axpy_prod (const matrix_expression<E1> &e1,
1.216 + const vector_expression<E2> &e2,
1.217 + V &v, packed_random_access_iterator_tag) {
1.218 + typedef typename E1::orientation_category orientation_category;
1.219 + return axpy_prod (e1, e2, v, packed_random_access_iterator_tag (), orientation_category ());
1.220 + }
1.221 +
1.222 +
1.223 + /** \brief computes <tt>v += A x</tt> or <tt>v = A x</tt> in an
1.224 + optimized fashion.
1.225 +
1.226 + \param e1 the matrix expression \c A
1.227 + \param e2 the vector expression \c x
1.228 + \param v the result vector \c v
1.229 + \param init a boolean parameter
1.230 +
1.231 + <tt>axpy_prod(A, x, v, init)</tt> implements the well known
1.232 + axpy-product. Setting \a init to \c true is equivalent to call
1.233 + <tt>v.clear()</tt> before <tt>axpy_prod</tt>. Currently \a init
1.234 + defaults to \c true, but this may change in the future.
1.235 +
1.236 + Up to now there are some specialisation for compressed
1.237 + matrices that give a large speed up compared to prod.
1.238 +
1.239 + \ingroup blas2
1.240 +
1.241 + \internal
1.242 +
1.243 + template parameters:
1.244 + \param V type of the result vector \c v
1.245 + \param E1 type of a matrix expression \c A
1.246 + \param E2 type of a vector expression \c x
1.247 + */
1.248 + template<class V, class E1, class E2>
1.249 + BOOST_UBLAS_INLINE
1.250 + V &
1.251 + axpy_prod (const matrix_expression<E1> &e1,
1.252 + const vector_expression<E2> &e2,
1.253 + V &v, bool init = true) {
1.254 + typedef typename V::value_type value_type;
1.255 + typedef typename E2::const_iterator::iterator_category iterator_category;
1.256 +
1.257 + if (init)
1.258 + v.assign (zero_vector<value_type> (e1 ().size1 ()));
1.259 +#if BOOST_UBLAS_TYPE_CHECK
1.260 + vector<value_type> cv (v);
1.261 + typedef typename type_traits<value_type>::real_type real_type;
1.262 + real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
1.263 + indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));
1.264 +#endif
1.265 + axpy_prod (e1, e2, v, iterator_category ());
1.266 +#if BOOST_UBLAS_TYPE_CHECK
1.267 + BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
1.268 +#endif
1.269 + return v;
1.270 + }
1.271 + template<class V, class E1, class E2>
1.272 + BOOST_UBLAS_INLINE
1.273 + V
1.274 + axpy_prod (const matrix_expression<E1> &e1,
1.275 + const vector_expression<E2> &e2) {
1.276 + typedef V vector_type;
1.277 +
1.278 + vector_type v (e1 ().size1 ());
1.279 + return axpy_prod (e1, e2, v, true);
1.280 + }
1.281 +
1.282 + template<class V, class E1, class T2, class IA2, class TA2>
1.283 + BOOST_UBLAS_INLINE
1.284 + V &
1.285 + axpy_prod (const vector_expression<E1> &e1,
1.286 + const compressed_matrix<T2, column_major, 0, IA2, TA2> &e2,
1.287 + V &v, column_major_tag) {
1.288 + typedef typename V::size_type size_type;
1.289 + typedef typename V::value_type value_type;
1.290 +
1.291 + for (size_type j = 0; j < e2.filled1 () -1; ++ j) {
1.292 + size_type begin = e2.index1_data () [j];
1.293 + size_type end = e2.index1_data () [j + 1];
1.294 + value_type t (v (j));
1.295 + for (size_type i = begin; i < end; ++ i)
1.296 + t += e2.value_data () [i] * e1 () (e2.index2_data () [i]);
1.297 + v (j) = t;
1.298 + }
1.299 + return v;
1.300 + }
1.301 +
1.302 + template<class V, class E1, class T2, class IA2, class TA2>
1.303 + BOOST_UBLAS_INLINE
1.304 + V &
1.305 + axpy_prod (const vector_expression<E1> &e1,
1.306 + const compressed_matrix<T2, row_major, 0, IA2, TA2> &e2,
1.307 + V &v, row_major_tag) {
1.308 + typedef typename V::size_type size_type;
1.309 +
1.310 + for (size_type i = 0; i < e2.filled1 () -1; ++ i) {
1.311 + size_type begin = e2.index1_data () [i];
1.312 + size_type end = e2.index1_data () [i + 1];
1.313 + for (size_type j = begin; j < end; ++ j)
1.314 + v (e2.index2_data () [j]) += e2.value_data () [j] * e1 () (i);
1.315 + }
1.316 + return v;
1.317 + }
1.318 +
1.319 + // Dispatcher
1.320 + template<class V, class E1, class T2, class L2, class IA2, class TA2>
1.321 + BOOST_UBLAS_INLINE
1.322 + V &
1.323 + axpy_prod (const vector_expression<E1> &e1,
1.324 + const compressed_matrix<T2, L2, 0, IA2, TA2> &e2,
1.325 + V &v, bool init = true) {
1.326 + typedef typename V::value_type value_type;
1.327 + typedef typename L2::orientation_category orientation_category;
1.328 +
1.329 + if (init)
1.330 + v.assign (zero_vector<value_type> (e2.size2 ()));
1.331 +#if BOOST_UBLAS_TYPE_CHECK
1.332 + vector<value_type> cv (v);
1.333 + typedef typename type_traits<value_type>::real_type real_type;
1.334 + real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
1.335 + indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));
1.336 +#endif
1.337 + axpy_prod (e1, e2, v, orientation_category ());
1.338 +#if BOOST_UBLAS_TYPE_CHECK
1.339 + BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
1.340 +#endif
1.341 + return v;
1.342 + }
1.343 + template<class V, class E1, class T2, class L2, class IA2, class TA2>
1.344 + BOOST_UBLAS_INLINE
1.345 + V
1.346 + axpy_prod (const vector_expression<E1> &e1,
1.347 + const compressed_matrix<T2, L2, 0, IA2, TA2> &e2) {
1.348 + typedef V vector_type;
1.349 +
1.350 + vector_type v (e2.size2 ());
1.351 + return axpy_prod (e1, e2, v, true);
1.352 + }
1.353 +
1.354 + template<class V, class E1, class E2>
1.355 + BOOST_UBLAS_INLINE
1.356 + V &
1.357 + axpy_prod (const vector_expression<E1> &e1,
1.358 + const matrix_expression<E2> &e2,
1.359 + V &v, packed_random_access_iterator_tag, column_major_tag) {
1.360 + typedef const E1 expression1_type;
1.361 + typedef const E2 expression2_type;
1.362 + typedef typename V::size_type size_type;
1.363 +
1.364 + typename expression2_type::const_iterator2 it2 (e2 ().begin2 ());
1.365 + typename expression2_type::const_iterator2 it2_end (e2 ().end2 ());
1.366 + while (it2 != it2_end) {
1.367 + size_type index2 (it2.index2 ());
1.368 +#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1.369 + typename expression2_type::const_iterator1 it1 (it2.begin ());
1.370 + typename expression2_type::const_iterator1 it1_end (it2.end ());
1.371 +#else
1.372 + typename expression2_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));
1.373 + typename expression2_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));
1.374 +#endif
1.375 + while (it1 != it1_end) {
1.376 + v (index2) += *it1 * e1 () (it1.index1 ());
1.377 + ++ it1;
1.378 + }
1.379 + ++ it2;
1.380 + }
1.381 + return v;
1.382 + }
1.383 +
1.384 + template<class V, class E1, class E2>
1.385 + BOOST_UBLAS_INLINE
1.386 + V &
1.387 + axpy_prod (const vector_expression<E1> &e1,
1.388 + const matrix_expression<E2> &e2,
1.389 + V &v, packed_random_access_iterator_tag, row_major_tag) {
1.390 + typedef const E1 expression1_type;
1.391 + typedef const E2 expression2_type;
1.392 + typedef typename V::size_type size_type;
1.393 +
1.394 + typename expression2_type::const_iterator1 it1 (e2 ().begin1 ());
1.395 + typename expression2_type::const_iterator1 it1_end (e2 ().end1 ());
1.396 + while (it1 != it1_end) {
1.397 + size_type index1 (it1.index1 ());
1.398 +#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1.399 + typename expression2_type::const_iterator2 it2 (it1.begin ());
1.400 + typename expression2_type::const_iterator2 it2_end (it1.end ());
1.401 +#else
1.402 + typename expression2_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));
1.403 + typename expression2_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));
1.404 +#endif
1.405 + while (it2 != it2_end) {
1.406 + v (it2.index2 ()) += *it2 * e1 () (index1);
1.407 + ++ it2;
1.408 + }
1.409 + ++ it1;
1.410 + }
1.411 + return v;
1.412 + }
1.413 +
1.414 + template<class V, class E1, class E2>
1.415 + BOOST_UBLAS_INLINE
1.416 + V &
1.417 + axpy_prod (const vector_expression<E1> &e1,
1.418 + const matrix_expression<E2> &e2,
1.419 + V &v, sparse_bidirectional_iterator_tag) {
1.420 + typedef const E1 expression1_type;
1.421 + typedef const E2 expression2_type;
1.422 + typedef typename V::size_type size_type;
1.423 +
1.424 + typename expression1_type::const_iterator it (e1 ().begin ());
1.425 + typename expression1_type::const_iterator it_end (e1 ().end ());
1.426 + while (it != it_end) {
1.427 + v.plus_assign (*it * row (e2 (), it.index ()));
1.428 + ++ it;
1.429 + }
1.430 + return v;
1.431 + }
1.432 +
1.433 + // Dispatcher
1.434 + template<class V, class E1, class E2>
1.435 + BOOST_UBLAS_INLINE
1.436 + V &
1.437 + axpy_prod (const vector_expression<E1> &e1,
1.438 + const matrix_expression<E2> &e2,
1.439 + V &v, packed_random_access_iterator_tag) {
1.440 + typedef typename E2::orientation_category orientation_category;
1.441 + return axpy_prod (e1, e2, v, packed_random_access_iterator_tag (), orientation_category ());
1.442 + }
1.443 +
1.444 +
1.445 + /** \brief computes <tt>v += A<sup>T</sup> x</tt> or <tt>v = A<sup>T</sup> x</tt> in an
1.446 + optimized fashion.
1.447 +
1.448 + \param e1 the vector expression \c x
1.449 + \param e2 the matrix expression \c A
1.450 + \param v the result vector \c v
1.451 + \param init a boolean parameter
1.452 +
1.453 + <tt>axpy_prod(x, A, v, init)</tt> implements the well known
1.454 + axpy-product. Setting \a init to \c true is equivalent to call
1.455 + <tt>v.clear()</tt> before <tt>axpy_prod</tt>. Currently \a init
1.456 + defaults to \c true, but this may change in the future.
1.457 +
1.458 + Up to now there are some specialisation for compressed
1.459 + matrices that give a large speed up compared to prod.
1.460 +
1.461 + \ingroup blas2
1.462 +
1.463 + \internal
1.464 +
1.465 + template parameters:
1.466 + \param V type of the result vector \c v
1.467 + \param E1 type of a vector expression \c x
1.468 + \param E2 type of a matrix expression \c A
1.469 + */
1.470 + template<class V, class E1, class E2>
1.471 + BOOST_UBLAS_INLINE
1.472 + V &
1.473 + axpy_prod (const vector_expression<E1> &e1,
1.474 + const matrix_expression<E2> &e2,
1.475 + V &v, bool init = true) {
1.476 + typedef typename V::value_type value_type;
1.477 + typedef typename E1::const_iterator::iterator_category iterator_category;
1.478 +
1.479 + if (init)
1.480 + v.assign (zero_vector<value_type> (e2 ().size2 ()));
1.481 +#if BOOST_UBLAS_TYPE_CHECK
1.482 + vector<value_type> cv (v);
1.483 + typedef typename type_traits<value_type>::real_type real_type;
1.484 + real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
1.485 + indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));
1.486 +#endif
1.487 + axpy_prod (e1, e2, v, iterator_category ());
1.488 +#if BOOST_UBLAS_TYPE_CHECK
1.489 + BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
1.490 +#endif
1.491 + return v;
1.492 + }
1.493 + template<class V, class E1, class E2>
1.494 + BOOST_UBLAS_INLINE
1.495 + V
1.496 + axpy_prod (const vector_expression<E1> &e1,
1.497 + const matrix_expression<E2> &e2) {
1.498 + typedef V vector_type;
1.499 +
1.500 + vector_type v (e2 ().size2 ());
1.501 + return axpy_prod (e1, e2, v, true);
1.502 + }
1.503 +
1.504 + template<class M, class E1, class E2, class TRI>
1.505 + BOOST_UBLAS_INLINE
1.506 + M &
1.507 + axpy_prod (const matrix_expression<E1> &e1,
1.508 + const matrix_expression<E2> &e2,
1.509 + M &m, TRI,
1.510 + dense_proxy_tag, row_major_tag) {
1.511 + typedef M matrix_type;
1.512 + typedef const E1 expression1_type;
1.513 + typedef const E2 expression2_type;
1.514 + typedef typename M::size_type size_type;
1.515 + typedef typename M::value_type value_type;
1.516 +
1.517 +#if BOOST_UBLAS_TYPE_CHECK
1.518 + matrix<value_type, row_major> cm (m);
1.519 + typedef typename type_traits<value_type>::real_type real_type;
1.520 + real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
1.521 + indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), row_major_tag ());
1.522 +#endif
1.523 + size_type size1 (e1 ().size1 ());
1.524 + size_type size2 (e1 ().size2 ());
1.525 + for (size_type i = 0; i < size1; ++ i)
1.526 + for (size_type j = 0; j < size2; ++ j)
1.527 + row (m, i).plus_assign (e1 () (i, j) * row (e2 (), j));
1.528 +#if BOOST_UBLAS_TYPE_CHECK
1.529 + BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
1.530 +#endif
1.531 + return m;
1.532 + }
1.533 + template<class M, class E1, class E2, class TRI>
1.534 + BOOST_UBLAS_INLINE
1.535 + M &
1.536 + axpy_prod (const matrix_expression<E1> &e1,
1.537 + const matrix_expression<E2> &e2,
1.538 + M &m, TRI,
1.539 + sparse_proxy_tag, row_major_tag) {
1.540 + typedef M matrix_type;
1.541 + typedef TRI triangular_restriction;
1.542 + typedef const E1 expression1_type;
1.543 + typedef const E2 expression2_type;
1.544 + typedef typename M::size_type size_type;
1.545 + typedef typename M::value_type value_type;
1.546 +
1.547 +#if BOOST_UBLAS_TYPE_CHECK
1.548 + matrix<value_type, row_major> cm (m);
1.549 + typedef typename type_traits<value_type>::real_type real_type;
1.550 + real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
1.551 + indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), row_major_tag ());
1.552 +#endif
1.553 + typename expression1_type::const_iterator1 it1 (e1 ().begin1 ());
1.554 + typename expression1_type::const_iterator1 it1_end (e1 ().end1 ());
1.555 + while (it1 != it1_end) {
1.556 +#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1.557 + typename expression1_type::const_iterator2 it2 (it1.begin ());
1.558 + typename expression1_type::const_iterator2 it2_end (it1.end ());
1.559 +#else
1.560 + typename expression1_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));
1.561 + typename expression1_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));
1.562 +#endif
1.563 + while (it2 != it2_end) {
1.564 + // row (m, it1.index1 ()).plus_assign (*it2 * row (e2 (), it2.index2 ()));
1.565 + matrix_row<expression2_type> mr (e2 (), it2.index2 ());
1.566 + typename matrix_row<expression2_type>::const_iterator itr (mr.begin ());
1.567 + typename matrix_row<expression2_type>::const_iterator itr_end (mr.end ());
1.568 + while (itr != itr_end) {
1.569 + if (triangular_restriction::other (it1.index1 (), itr.index ()))
1.570 + m (it1.index1 (), itr.index ()) += *it2 * *itr;
1.571 + ++ itr;
1.572 + }
1.573 + ++ it2;
1.574 + }
1.575 + ++ it1;
1.576 + }
1.577 +#if BOOST_UBLAS_TYPE_CHECK
1.578 + BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
1.579 +#endif
1.580 + return m;
1.581 + }
1.582 +
1.583 + template<class M, class E1, class E2, class TRI>
1.584 + BOOST_UBLAS_INLINE
1.585 + M &
1.586 + axpy_prod (const matrix_expression<E1> &e1,
1.587 + const matrix_expression<E2> &e2,
1.588 + M &m, TRI,
1.589 + dense_proxy_tag, column_major_tag) {
1.590 + typedef M matrix_type;
1.591 + typedef const E1 expression1_type;
1.592 + typedef const E2 expression2_type;
1.593 + typedef typename M::size_type size_type;
1.594 + typedef typename M::value_type value_type;
1.595 +
1.596 +#if BOOST_UBLAS_TYPE_CHECK
1.597 + matrix<value_type, column_major> cm (m);
1.598 + typedef typename type_traits<value_type>::real_type real_type;
1.599 + real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
1.600 + indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), column_major_tag ());
1.601 +#endif
1.602 + size_type size1 (e2 ().size1 ());
1.603 + size_type size2 (e2 ().size2 ());
1.604 + for (size_type j = 0; j < size2; ++ j)
1.605 + for (size_type i = 0; i < size1; ++ i)
1.606 + column (m, j).plus_assign (e2 () (i, j) * column (e1 (), i));
1.607 +#if BOOST_UBLAS_TYPE_CHECK
1.608 + BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
1.609 +#endif
1.610 + return m;
1.611 + }
1.612 + template<class M, class E1, class E2, class TRI>
1.613 + BOOST_UBLAS_INLINE
1.614 + M &
1.615 + axpy_prod (const matrix_expression<E1> &e1,
1.616 + const matrix_expression<E2> &e2,
1.617 + M &m, TRI,
1.618 + sparse_proxy_tag, column_major_tag) {
1.619 + typedef M matrix_type;
1.620 + typedef TRI triangular_restriction;
1.621 + typedef const E1 expression1_type;
1.622 + typedef const E2 expression2_type;
1.623 + typedef typename M::size_type size_type;
1.624 + typedef typename M::value_type value_type;
1.625 +
1.626 +#if BOOST_UBLAS_TYPE_CHECK
1.627 + matrix<value_type, column_major> cm (m);
1.628 + typedef typename type_traits<value_type>::real_type real_type;
1.629 + real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
1.630 + indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), column_major_tag ());
1.631 +#endif
1.632 + typename expression2_type::const_iterator2 it2 (e2 ().begin2 ());
1.633 + typename expression2_type::const_iterator2 it2_end (e2 ().end2 ());
1.634 + while (it2 != it2_end) {
1.635 +#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1.636 + typename expression2_type::const_iterator1 it1 (it2.begin ());
1.637 + typename expression2_type::const_iterator1 it1_end (it2.end ());
1.638 +#else
1.639 + typename expression2_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));
1.640 + typename expression2_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));
1.641 +#endif
1.642 + while (it1 != it1_end) {
1.643 + // column (m, it2.index2 ()).plus_assign (*it1 * column (e1 (), it1.index1 ()));
1.644 + matrix_column<expression1_type> mc (e1 (), it1.index1 ());
1.645 + typename matrix_column<expression1_type>::const_iterator itc (mc.begin ());
1.646 + typename matrix_column<expression1_type>::const_iterator itc_end (mc.end ());
1.647 + while (itc != itc_end) {
1.648 + if (triangular_restriction::functor_type ().other (itc.index (), it2.index2 ()))
1.649 + m (itc.index (), it2.index2 ()) += *it1 * *itc;
1.650 + ++ itc;
1.651 + }
1.652 + ++ it1;
1.653 + }
1.654 + ++ it2;
1.655 + }
1.656 +#if BOOST_UBLAS_TYPE_CHECK
1.657 + BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
1.658 +#endif
1.659 + return m;
1.660 + }
1.661 +
1.662 + // Dispatcher
1.663 + template<class M, class E1, class E2, class TRI>
1.664 + BOOST_UBLAS_INLINE
1.665 + M &
1.666 + axpy_prod (const matrix_expression<E1> &e1,
1.667 + const matrix_expression<E2> &e2,
1.668 + M &m, TRI, bool init = true) {
1.669 + typedef typename M::value_type value_type;
1.670 + typedef typename M::storage_category storage_category;
1.671 + typedef typename M::orientation_category orientation_category;
1.672 + typedef TRI triangular_restriction;
1.673 +
1.674 + if (init)
1.675 + m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
1.676 + return axpy_prod (e1, e2, m, triangular_restriction (), storage_category (), orientation_category ());
1.677 + }
1.678 + template<class M, class E1, class E2, class TRI>
1.679 + BOOST_UBLAS_INLINE
1.680 + M
1.681 + axpy_prod (const matrix_expression<E1> &e1,
1.682 + const matrix_expression<E2> &e2,
1.683 + TRI) {
1.684 + typedef M matrix_type;
1.685 + typedef TRI triangular_restriction;
1.686 +
1.687 + matrix_type m (e1 ().size1 (), e2 ().size2 ());
1.688 + return axpy_prod (e1, e2, m, triangular_restriction (), true);
1.689 + }
1.690 +
1.691 + /** \brief computes <tt>M += A X</tt> or <tt>M = A X</tt> in an
1.692 + optimized fashion.
1.693 +
1.694 + \param e1 the matrix expression \c A
1.695 + \param e2 the matrix expression \c X
1.696 + \param m the result matrix \c M
1.697 + \param init a boolean parameter
1.698 +
1.699 + <tt>axpy_prod(A, X, M, init)</tt> implements the well known
1.700 + axpy-product. Setting \a init to \c true is equivalent to call
1.701 + <tt>M.clear()</tt> before <tt>axpy_prod</tt>. Currently \a init
1.702 + defaults to \c true, but this may change in the future.
1.703 +
1.704 + Up to now there are no specialisations.
1.705 +
1.706 + \ingroup blas3
1.707 +
1.708 + \internal
1.709 +
1.710 + template parameters:
1.711 + \param M type of the result matrix \c M
1.712 + \param E1 type of a matrix expression \c A
1.713 + \param E2 type of a matrix expression \c X
1.714 + */
1.715 + template<class M, class E1, class E2>
1.716 + BOOST_UBLAS_INLINE
1.717 + M &
1.718 + axpy_prod (const matrix_expression<E1> &e1,
1.719 + const matrix_expression<E2> &e2,
1.720 + M &m, bool init = true) {
1.721 + typedef typename M::value_type value_type;
1.722 + typedef typename M::storage_category storage_category;
1.723 + typedef typename M::orientation_category orientation_category;
1.724 +
1.725 + if (init)
1.726 + m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
1.727 + return axpy_prod (e1, e2, m, full (), storage_category (), orientation_category ());
1.728 + }
1.729 + template<class M, class E1, class E2>
1.730 + BOOST_UBLAS_INLINE
1.731 + M
1.732 + axpy_prod (const matrix_expression<E1> &e1,
1.733 + const matrix_expression<E2> &e2) {
1.734 + typedef M matrix_type;
1.735 +
1.736 + matrix_type m (e1 ().size1 (), e2 ().size2 ());
1.737 + return axpy_prod (e1, e2, m, full (), true);
1.738 + }
1.739 +
1.740 +
1.741 + template<class M, class E1, class E2>
1.742 + BOOST_UBLAS_INLINE
1.743 + M &
1.744 + opb_prod (const matrix_expression<E1> &e1,
1.745 + const matrix_expression<E2> &e2,
1.746 + M &m,
1.747 + dense_proxy_tag, row_major_tag) {
1.748 + typedef M matrix_type;
1.749 + typedef const E1 expression1_type;
1.750 + typedef const E2 expression2_type;
1.751 + typedef typename M::size_type size_type;
1.752 + typedef typename M::value_type value_type;
1.753 +
1.754 +#if BOOST_UBLAS_TYPE_CHECK
1.755 + matrix<value_type, row_major> cm (m);
1.756 + typedef typename type_traits<value_type>::real_type real_type;
1.757 + real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
1.758 + indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), row_major_tag ());
1.759 +#endif
1.760 + size_type size (BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ()));
1.761 + for (size_type k = 0; k < size; ++ k) {
1.762 + vector<value_type> ce1 (column (e1 (), k));
1.763 + vector<value_type> re2 (row (e2 (), k));
1.764 + m.plus_assign (outer_prod (ce1, re2));
1.765 + }
1.766 +#if BOOST_UBLAS_TYPE_CHECK
1.767 + BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
1.768 +#endif
1.769 + return m;
1.770 + }
1.771 +
1.772 + template<class M, class E1, class E2>
1.773 + BOOST_UBLAS_INLINE
1.774 + M &
1.775 + opb_prod (const matrix_expression<E1> &e1,
1.776 + const matrix_expression<E2> &e2,
1.777 + M &m,
1.778 + dense_proxy_tag, column_major_tag) {
1.779 + typedef M matrix_type;
1.780 + typedef const E1 expression1_type;
1.781 + typedef const E2 expression2_type;
1.782 + typedef typename M::size_type size_type;
1.783 + typedef typename M::value_type value_type;
1.784 +
1.785 +#if BOOST_UBLAS_TYPE_CHECK
1.786 + matrix<value_type, column_major> cm (m);
1.787 + typedef typename type_traits<value_type>::real_type real_type;
1.788 + real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
1.789 + indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), column_major_tag ());
1.790 +#endif
1.791 + size_type size (BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ()));
1.792 + for (size_type k = 0; k < size; ++ k) {
1.793 + vector<value_type> ce1 (column (e1 (), k));
1.794 + vector<value_type> re2 (row (e2 (), k));
1.795 + m.plus_assign (outer_prod (ce1, re2));
1.796 + }
1.797 +#if BOOST_UBLAS_TYPE_CHECK
1.798 + BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
1.799 +#endif
1.800 + return m;
1.801 + }
1.802 +
1.803 + // Dispatcher
1.804 +
1.805 + /** \brief computes <tt>M += A X</tt> or <tt>M = A X</tt> in an
1.806 + optimized fashion.
1.807 +
1.808 + \param e1 the matrix expression \c A
1.809 + \param e2 the matrix expression \c X
1.810 + \param m the result matrix \c M
1.811 + \param init a boolean parameter
1.812 +
1.813 + <tt>opb_prod(A, X, M, init)</tt> implements the well known
1.814 + axpy-product. Setting \a init to \c true is equivalent to call
1.815 + <tt>M.clear()</tt> before <tt>opb_prod</tt>. Currently \a init
1.816 + defaults to \c true, but this may change in the future.
1.817 +
1.818 + This function may give a speedup if \c A has less columns than
1.819 + rows, because the product is computed as a sum of outer
1.820 + products.
1.821 +
1.822 + \ingroup blas3
1.823 +
1.824 + \internal
1.825 +
1.826 + template parameters:
1.827 + \param M type of the result matrix \c M
1.828 + \param E1 type of a matrix expression \c A
1.829 + \param E2 type of a matrix expression \c X
1.830 + */
1.831 + template<class M, class E1, class E2>
1.832 + BOOST_UBLAS_INLINE
1.833 + M &
1.834 + opb_prod (const matrix_expression<E1> &e1,
1.835 + const matrix_expression<E2> &e2,
1.836 + M &m, bool init = true) {
1.837 + typedef typename M::value_type value_type;
1.838 + typedef typename M::storage_category storage_category;
1.839 + typedef typename M::orientation_category orientation_category;
1.840 +
1.841 + if (init)
1.842 + m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
1.843 + return opb_prod (e1, e2, m, storage_category (), orientation_category ());
1.844 + }
1.845 + template<class M, class E1, class E2>
1.846 + BOOST_UBLAS_INLINE
1.847 + M
1.848 + opb_prod (const matrix_expression<E1> &e1,
1.849 + const matrix_expression<E2> &e2) {
1.850 + typedef M matrix_type;
1.851 +
1.852 + matrix_type m (e1 ().size1 (), e2 ().size2 ());
1.853 + return opb_prod (e1, e2, m, true);
1.854 + }
1.855 +
1.856 +}}}
1.857 +
1.858 +#endif