numeric/ublas/operation.hpp

0001 //
0002 //  Copyright (c) 2000-2002
0003 //  Joerg Walter, Mathias Koch
0004 //
0005 //  Distributed under the Boost Software License, Version 1.0. (See
0006 //  accompanying file LICENSE_1_0.txt or copy at
0007 //  http://www.boost.org/LICENSE_1_0.txt)
0008 //
0009 //  The authors gratefully acknowledge the support of
0010 //  GeNeSys mbH & Co. KG in producing this work.
0011 //
0012
0013 #ifndef _BOOST_UBLAS_OPERATION_
0014 #define _BOOST_UBLAS_OPERATION_
0015
0016 #include <boost/numeric/ublas/matrix_proxy.hpp>
0017
0018 /** \file operation.hpp
0019  *  \brief This file contains some specialized products.
0020  */
0021
0022 // axpy-based products
0023 // Alexei Novakov had a lot of ideas to improve these. Thanks.
0024 // Hendrik Kueck proposed some new kernel. Thanks again.
0025
0026 namespace boost { namespace numeric { namespace ublas {
0027
0028     template<class V, class T1, class L1, class IA1, class TA1, class E2>
0029     BOOST_UBLAS_INLINE
0030     V &
0031     axpy_prod (const compressed_matrix<T1, L1, 0, IA1, TA1> &e1,
0032                const vector_expression<E2> &e2,
0033                V &v, row_major_tag) {
0034         typedef typename V::size_type size_type;
0035         typedef typename V::value_type value_type;
0036
0037         for (size_type i = 0; i < e1.filled1 () -1; ++ i) {
0038             size_type begin = e1.index1_data () [i];
0039             size_type end = e1.index1_data () [i + 1];
0040             value_type t (v (i));
0041             for (size_type j = begin; j < end; ++ j)
0042                 t += e1.value_data () [j] * e2 () (e1.index2_data () [j]);
0043             v (i) = t;
0044         }
0045         return v;
0046     }
0047
0048     template<class V, class T1, class L1, class IA1, class TA1, class E2>
0049     BOOST_UBLAS_INLINE
0050     V &
0051     axpy_prod (const compressed_matrix<T1, L1, 0, IA1, TA1> &e1,
0052                const vector_expression<E2> &e2,
0053                V &v, column_major_tag) {
0054         typedef typename V::size_type size_type;
0055
0056         for (size_type j = 0; j < e1.filled1 () -1; ++ j) {
0057             size_type begin = e1.index1_data () [j];
0058             size_type end = e1.index1_data () [j + 1];
0059             for (size_type i = begin; i < end; ++ i)
0060                 v (e1.index2_data () [i]) += e1.value_data () [i] * e2 () (j);
0061         }
0062         return v;
0063     }
0064
0065     // Dispatcher
0066     template<class V, class T1, class L1, class IA1, class TA1, class E2>
0067     BOOST_UBLAS_INLINE
0068     V &
0069     axpy_prod (const compressed_matrix<T1, L1, 0, IA1, TA1> &e1,
0070                const vector_expression<E2> &e2,
0071                V &v, bool init = true) {
0072         typedef typename V::value_type value_type;
0073         typedef typename L1::orientation_category orientation_category;
0074
0075         if (init)
0076             v.assign (zero_vector<value_type> (e1.size1 ()));
0077 #if BOOST_UBLAS_TYPE_CHECK
0078         vector<value_type> cv (v);
0079         typedef typename type_traits<value_type>::real_type real_type;
0080         real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
0081         indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));
0082 #endif
0083         axpy_prod (e1, e2, v, orientation_category ());
0084 #if BOOST_UBLAS_TYPE_CHECK
0085         BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
0086 #endif
0087         return v;
0088     }
0089     template<class V, class T1, class L1, class IA1, class TA1, class E2>
0090     BOOST_UBLAS_INLINE
0091     V
0092     axpy_prod (const compressed_matrix<T1, L1, 0, IA1, TA1> &e1,
0093                const vector_expression<E2> &e2) {
0094         typedef V vector_type;
0095
0096         vector_type v (e1.size1 ());
0097         return axpy_prod (e1, e2, v, true);
0098     }
0099
0100     template<class V, class T1, class L1, class IA1, class TA1, class E2>
0101     BOOST_UBLAS_INLINE
0102     V &
0103     axpy_prod (const coordinate_matrix<T1, L1, 0, IA1, TA1> &e1,
0104                const vector_expression<E2> &e2,
0105                V &v, bool init = true) {
0106         typedef typename V::size_type size_type;
0107         typedef typename V::value_type value_type;
0108         typedef L1 layout_type;
0109
0110         size_type size1 = e1.size1();
0111         size_type size2 = e1.size2();
0112
0113         if (init) {
0114             noalias(v) = zero_vector<value_type>(size1);
0115         }
0116
0117         for (size_type i = 0; i < e1.nnz(); ++i) {
0118             size_type row_index = layout_type::index_M( e1.index1_data () [i], e1.index2_data () [i] );
0119             size_type col_index = layout_type::index_m( e1.index1_data () [i], e1.index2_data () [i] );
0120             v( row_index ) += e1.value_data () [i] * e2 () (col_index);
0121         }
0122         return v;
0123     }
0124
0125     template<class V, class E1, class E2>
0126     BOOST_UBLAS_INLINE
0127     V &
0128     axpy_prod (const matrix_expression<E1> &e1,
0129                const vector_expression<E2> &e2,
0130                V &v, packed_random_access_iterator_tag, row_major_tag) {
0131         typedef const E1 expression1_type;
0132         typedef typename V::size_type size_type;
0133
0134         typename expression1_type::const_iterator1 it1 (e1 ().begin1 ());
0135         typename expression1_type::const_iterator1 it1_end (e1 ().end1 ());
0136         while (it1 != it1_end) {
0137             size_type index1 (it1.index1 ());
0138 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
0139             typename expression1_type::const_iterator2 it2 (it1.begin ());
0140             typename expression1_type::const_iterator2 it2_end (it1.end ());
0141 #else
0142             typename expression1_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));
0143             typename expression1_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));
0144 #endif
0145             while (it2 != it2_end) {
0146                 v (index1) += *it2 * e2 () (it2.index2 ());
0147                 ++ it2;
0148             }
0149             ++ it1;
0150         }
0151         return v;
0152     }
0153
0154     template<class V, class E1, class E2>
0155     BOOST_UBLAS_INLINE
0156     V &
0157     axpy_prod (const matrix_expression<E1> &e1,
0158                const vector_expression<E2> &e2,
0159                V &v, packed_random_access_iterator_tag, column_major_tag) {
0160         typedef const E1 expression1_type;
0161         typedef typename V::size_type size_type;
0162
0163         typename expression1_type::const_iterator2 it2 (e1 ().begin2 ());
0164         typename expression1_type::const_iterator2 it2_end (e1 ().end2 ());
0165         while (it2 != it2_end) {
0166             size_type index2 (it2.index2 ());
0167 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
0168             typename expression1_type::const_iterator1 it1 (it2.begin ());
0169             typename expression1_type::const_iterator1 it1_end (it2.end ());
0170 #else
0171             typename expression1_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));
0172             typename expression1_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));
0173 #endif
0174             while (it1 != it1_end) {
0175                 v (it1.index1 ()) += *it1 * e2 () (index2);
0176                 ++ it1;
0177             }
0178             ++ it2;
0179         }
0180         return v;
0181     }
0182
0183     template<class V, class E1, class E2>
0184     BOOST_UBLAS_INLINE
0185     V &
0186     axpy_prod (const matrix_expression<E1> &e1,
0187                const vector_expression<E2> &e2,
0188                V &v, sparse_bidirectional_iterator_tag) {
0189         typedef const E2 expression2_type;
0190
0191         typename expression2_type::const_iterator it (e2 ().begin ());
0192         typename expression2_type::const_iterator it_end (e2 ().end ());
0193         while (it != it_end) {
0194             v.plus_assign (column (e1 (), it.index ()) * *it);
0195             ++ it;
0196         }
0197         return v;
0198     }
0199
0200     // Dispatcher
0201     template<class V, class E1, class E2>
0202     BOOST_UBLAS_INLINE
0203     V &
0204     axpy_prod (const matrix_expression<E1> &e1,
0205                const vector_expression<E2> &e2,
0206                V &v, packed_random_access_iterator_tag) {
0207         typedef typename E1::orientation_category orientation_category;
0208         return axpy_prod (e1, e2, v, packed_random_access_iterator_tag (), orientation_category ());
0209     }
0210
0211
0212   /** \brief computes <tt>v += A x</tt> or <tt>v = A x</tt> in an
0213           optimized fashion.
0214
0215           \param e1 the matrix expression \c A
0216           \param e2 the vector expression \c x
0217           \param v  the result vector \c v
0218           \param init a boolean parameter
0219
0220           <tt>axpy_prod(A, x, v, init)</tt> implements the well known
0221           axpy-product.  Setting \a init to \c true is equivalent to call
0222           <tt>v.clear()</tt> before <tt>axpy_prod</tt>. Currently \a init
0223           defaults to \c true, but this may change in the future.
0224
0225           Up to now there are some specialisation for compressed
0226           matrices that give a large speed up compared to prod.
0227
0228           \ingroup blas2
0229
0230           \internal
0231
0232           template parameters:
0233           \param V type of the result vector \c v
0234           \param E1 type of a matrix expression \c A
0235           \param E2 type of a vector expression \c x
0236   */
0237     template<class V, class E1, class E2>
0238     BOOST_UBLAS_INLINE
0239     V &
0240     axpy_prod (const matrix_expression<E1> &e1,
0241                const vector_expression<E2> &e2,
0242                V &v, bool init = true) {
0243         typedef typename V::value_type value_type;
0244         typedef typename E2::const_iterator::iterator_category iterator_category;
0245
0246         if (init)
0247             v.assign (zero_vector<value_type> (e1 ().size1 ()));
0248 #if BOOST_UBLAS_TYPE_CHECK
0249         vector<value_type> cv (v);
0250         typedef typename type_traits<value_type>::real_type real_type;
0251         real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
0252         indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));
0253 #endif
0254         axpy_prod (e1, e2, v, iterator_category ());
0255 #if BOOST_UBLAS_TYPE_CHECK
0256         BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
0257 #endif
0258         return v;
0259     }
0260     template<class V, class E1, class E2>
0261     BOOST_UBLAS_INLINE
0262     V
0263     axpy_prod (const matrix_expression<E1> &e1,
0264                const vector_expression<E2> &e2) {
0265         typedef V vector_type;
0266
0267         vector_type v (e1 ().size1 ());
0268         return axpy_prod (e1, e2, v, true);
0269     }
0270
0271     template<class V, class E1, class T2, class IA2, class TA2>
0272     BOOST_UBLAS_INLINE
0273     V &
0274     axpy_prod (const vector_expression<E1> &e1,
0275                const compressed_matrix<T2, column_major, 0, IA2, TA2> &e2,
0276                V &v, column_major_tag) {
0277         typedef typename V::size_type size_type;
0278         typedef typename V::value_type value_type;
0279
0280         for (size_type j = 0; j < e2.filled1 () -1; ++ j) {
0281             size_type begin = e2.index1_data () [j];
0282             size_type end = e2.index1_data () [j + 1];
0283             value_type t (v (j));
0284             for (size_type i = begin; i < end; ++ i)
0285                 t += e2.value_data () [i] * e1 () (e2.index2_data () [i]);
0286             v (j) = t;
0287         }
0288         return v;
0289     }
0290
0291     template<class V, class E1, class T2, class IA2, class TA2>
0292     BOOST_UBLAS_INLINE
0293     V &
0294     axpy_prod (const vector_expression<E1> &e1,
0295                const compressed_matrix<T2, row_major, 0, IA2, TA2> &e2,
0296                V &v, row_major_tag) {
0297         typedef typename V::size_type size_type;
0298
0299         for (size_type i = 0; i < e2.filled1 () -1; ++ i) {
0300             size_type begin = e2.index1_data () [i];
0301             size_type end = e2.index1_data () [i + 1];
0302             for (size_type j = begin; j < end; ++ j)
0303                 v (e2.index2_data () [j]) += e2.value_data () [j] * e1 () (i);
0304         }
0305         return v;
0306     }
0307
0308     // Dispatcher
0309     template<class V, class E1, class T2, class L2, class IA2, class TA2>
0310     BOOST_UBLAS_INLINE
0311     V &
0312     axpy_prod (const vector_expression<E1> &e1,
0313                const compressed_matrix<T2, L2, 0, IA2, TA2> &e2,
0314                V &v, bool init = true) {
0315         typedef typename V::value_type value_type;
0316         typedef typename L2::orientation_category orientation_category;
0317
0318         if (init)
0319             v.assign (zero_vector<value_type> (e2.size2 ()));
0320 #if BOOST_UBLAS_TYPE_CHECK
0321         vector<value_type> cv (v);
0322         typedef typename type_traits<value_type>::real_type real_type;
0323         real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
0324         indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));
0325 #endif
0326         axpy_prod (e1, e2, v, orientation_category ());
0327 #if BOOST_UBLAS_TYPE_CHECK
0328         BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
0329 #endif
0330         return v;
0331     }
0332     template<class V, class E1, class T2, class L2, class IA2, class TA2>
0333     BOOST_UBLAS_INLINE
0334     V
0335     axpy_prod (const vector_expression<E1> &e1,
0336                const compressed_matrix<T2, L2, 0, IA2, TA2> &e2) {
0337         typedef V vector_type;
0338
0339         vector_type v (e2.size2 ());
0340         return axpy_prod (e1, e2, v, true);
0341     }
0342
0343     template<class V, class E1, class E2>
0344     BOOST_UBLAS_INLINE
0345     V &
0346     axpy_prod (const vector_expression<E1> &e1,
0347                const matrix_expression<E2> &e2,
0348                V &v, packed_random_access_iterator_tag, column_major_tag) {
0349         typedef const E2 expression2_type;
0350         typedef typename V::size_type size_type;
0351
0352         typename expression2_type::const_iterator2 it2 (e2 ().begin2 ());
0353         typename expression2_type::const_iterator2 it2_end (e2 ().end2 ());
0354         while (it2 != it2_end) {
0355             size_type index2 (it2.index2 ());
0356 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
0357             typename expression2_type::const_iterator1 it1 (it2.begin ());
0358             typename expression2_type::const_iterator1 it1_end (it2.end ());
0359 #else
0360             typename expression2_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));
0361             typename expression2_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));
0362 #endif
0363             while (it1 != it1_end) {
0364                 v (index2) += *it1 * e1 () (it1.index1 ());
0365                 ++ it1;
0366             }
0367             ++ it2;
0368         }
0369         return v;
0370     }
0371
0372     template<class V, class E1, class E2>
0373     BOOST_UBLAS_INLINE
0374     V &
0375     axpy_prod (const vector_expression<E1> &e1,
0376                const matrix_expression<E2> &e2,
0377                V &v, packed_random_access_iterator_tag, row_major_tag) {
0378         typedef const E2 expression2_type;
0379         typedef typename V::size_type size_type;
0380
0381         typename expression2_type::const_iterator1 it1 (e2 ().begin1 ());
0382         typename expression2_type::const_iterator1 it1_end (e2 ().end1 ());
0383         while (it1 != it1_end) {
0384             size_type index1 (it1.index1 ());
0385 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
0386             typename expression2_type::const_iterator2 it2 (it1.begin ());
0387             typename expression2_type::const_iterator2 it2_end (it1.end ());
0388 #else
0389             typename expression2_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));
0390             typename expression2_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));
0391 #endif
0392             while (it2 != it2_end) {
0393                 v (it2.index2 ()) += *it2 * e1 () (index1);
0394                 ++ it2;
0395             }
0396             ++ it1;
0397         }
0398         return v;
0399     }
0400
0401     template<class V, class E1, class E2>
0402     BOOST_UBLAS_INLINE
0403     V &
0404     axpy_prod (const vector_expression<E1> &e1,
0405                const matrix_expression<E2> &e2,
0406                V &v, sparse_bidirectional_iterator_tag) {
0407         typedef const E1 expression1_type;
0408
0409         typename expression1_type::const_iterator it (e1 ().begin ());
0410         typename expression1_type::const_iterator it_end (e1 ().end ());
0411         while (it != it_end) {
0412             v.plus_assign (*it * row (e2 (), it.index ()));
0413             ++ it;
0414         }
0415         return v;
0416     }
0417
0418     // Dispatcher
0419     template<class V, class E1, class E2>
0420     BOOST_UBLAS_INLINE
0421     V &
0422     axpy_prod (const vector_expression<E1> &e1,
0423                const matrix_expression<E2> &e2,
0424                V &v, packed_random_access_iterator_tag) {
0425         typedef typename E2::orientation_category orientation_category;
0426         return axpy_prod (e1, e2, v, packed_random_access_iterator_tag (), orientation_category ());
0427     }
0428
0429
0430   /** \brief computes <tt>v += A<sup>T</sup> x</tt> or <tt>v = A<sup>T</sup> x</tt> in an
0431           optimized fashion.
0432
0433           \param e1 the vector expression \c x
0434           \param e2 the matrix expression \c A
0435           \param v  the result vector \c v
0436           \param init a boolean parameter
0437
0438           <tt>axpy_prod(x, A, v, init)</tt> implements the well known
0439           axpy-product.  Setting \a init to \c true is equivalent to call
0440           <tt>v.clear()</tt> before <tt>axpy_prod</tt>. Currently \a init
0441           defaults to \c true, but this may change in the future.
0442
0443           Up to now there are some specialisation for compressed
0444           matrices that give a large speed up compared to prod.
0445
0446           \ingroup blas2
0447
0448           \internal
0449
0450           template parameters:
0451           \param V type of the result vector \c v
0452           \param E1 type of a vector expression \c x
0453           \param E2 type of a matrix expression \c A
0454   */
0455     template<class V, class E1, class E2>
0456     BOOST_UBLAS_INLINE
0457     V &
0458     axpy_prod (const vector_expression<E1> &e1,
0459                const matrix_expression<E2> &e2,
0460                V &v, bool init = true) {
0461         typedef typename V::value_type value_type;
0462         typedef typename E1::const_iterator::iterator_category iterator_category;
0463
0464         if (init)
0465             v.assign (zero_vector<value_type> (e2 ().size2 ()));
0466 #if BOOST_UBLAS_TYPE_CHECK
0467         vector<value_type> cv (v);
0468         typedef typename type_traits<value_type>::real_type real_type;
0469         real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
0470         indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));
0471 #endif
0472         axpy_prod (e1, e2, v, iterator_category ());
0473 #if BOOST_UBLAS_TYPE_CHECK
0474         BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
0475 #endif
0476         return v;
0477     }
0478     template<class V, class E1, class E2>
0479     BOOST_UBLAS_INLINE
0480     V
0481     axpy_prod (const vector_expression<E1> &e1,
0482                const matrix_expression<E2> &e2) {
0483         typedef V vector_type;
0484
0485         vector_type v (e2 ().size2 ());
0486         return axpy_prod (e1, e2, v, true);
0487     }
0488
0489     template<class M, class E1, class E2, class TRI>
0490     BOOST_UBLAS_INLINE
0491     M &
0492     axpy_prod (const matrix_expression<E1> &e1,
0493                const matrix_expression<E2> &e2,
0494                M &m, TRI,
0495                dense_proxy_tag, row_major_tag) {
0496
0497         typedef typename M::size_type size_type;
0498
0499 #if BOOST_UBLAS_TYPE_CHECK
0500         typedef typename M::value_type value_type;
0501         matrix<value_type, row_major> cm (m);
0502         typedef typename type_traits<value_type>::real_type real_type;
0503         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
0504         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), row_major_tag ());
0505 #endif
0506         size_type size1 (e1 ().size1 ());
0507         size_type size2 (e1 ().size2 ());
0508         for (size_type i = 0; i < size1; ++ i)
0509             for (size_type j = 0; j < size2; ++ j)
0510                 row (m, i).plus_assign (e1 () (i, j) * row (e2 (), j));
0511 #if BOOST_UBLAS_TYPE_CHECK
0512         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
0513 #endif
0514         return m;
0515     }
0516     template<class M, class E1, class E2, class TRI>
0517     BOOST_UBLAS_INLINE
0518     M &
0519     axpy_prod (const matrix_expression<E1> &e1,
0520                const matrix_expression<E2> &e2,
0521                M &m, TRI,
0522                sparse_proxy_tag, row_major_tag) {
0523
0524         typedef TRI triangular_restriction;
0525         typedef const E1 expression1_type;
0526         typedef const E2 expression2_type;
0527
0528 #if BOOST_UBLAS_TYPE_CHECK
0529         typedef typename M::value_type value_type;
0530         matrix<value_type, row_major> cm (m);
0531         typedef typename type_traits<value_type>::real_type real_type;
0532         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
0533         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), row_major_tag ());
0534 #endif
0535         typename expression1_type::const_iterator1 it1 (e1 ().begin1 ());
0536         typename expression1_type::const_iterator1 it1_end (e1 ().end1 ());
0537         while (it1 != it1_end) {
0538 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
0539             typename expression1_type::const_iterator2 it2 (it1.begin ());
0540             typename expression1_type::const_iterator2 it2_end (it1.end ());
0541 #else
0542             typename expression1_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));
0543             typename expression1_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));
0544 #endif
0545             while (it2 != it2_end) {
0546                 // row (m, it1.index1 ()).plus_assign (*it2 * row (e2 (), it2.index2 ()));
0547                 matrix_row<expression2_type> mr (e2 (), it2.index2 ());
0548                 typename matrix_row<expression2_type>::const_iterator itr (mr.begin ());
0549                 typename matrix_row<expression2_type>::const_iterator itr_end (mr.end ());
0550                 while (itr != itr_end) {
0551                     if (triangular_restriction::other (it1.index1 (), itr.index ()))
0552                         m (it1.index1 (), itr.index ()) += *it2 * *itr;
0553                     ++ itr;
0554                 }
0555                 ++ it2;
0556             }
0557             ++ it1;
0558         }
0559 #if BOOST_UBLAS_TYPE_CHECK
0560         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
0561 #endif
0562         return m;
0563     }
0564
0565     template<class M, class E1, class E2, class TRI>
0566     BOOST_UBLAS_INLINE
0567     M &
0568     axpy_prod (const matrix_expression<E1> &e1,
0569                const matrix_expression<E2> &e2,
0570                M &m, TRI,
0571                dense_proxy_tag, column_major_tag) {
0572         typedef typename M::size_type size_type;
0573
0574 #if BOOST_UBLAS_TYPE_CHECK
0575         typedef typename M::value_type value_type;
0576         matrix<value_type, column_major> cm (m);
0577         typedef typename type_traits<value_type>::real_type real_type;
0578         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
0579         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), column_major_tag ());
0580 #endif
0581         size_type size1 (e2 ().size1 ());
0582         size_type size2 (e2 ().size2 ());
0583         for (size_type j = 0; j < size2; ++ j)
0584             for (size_type i = 0; i < size1; ++ i)
0585                 column (m, j).plus_assign (e2 () (i, j) * column (e1 (), i));
0586 #if BOOST_UBLAS_TYPE_CHECK
0587         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
0588 #endif
0589         return m;
0590     }
0591     template<class M, class E1, class E2, class TRI>
0592     BOOST_UBLAS_INLINE
0593     M &
0594     axpy_prod (const matrix_expression<E1> &e1,
0595                const matrix_expression<E2> &e2,
0596                M &m, TRI,
0597                sparse_proxy_tag, column_major_tag) {
0598         typedef TRI triangular_restriction;
0599         typedef const E1 expression1_type;
0600         typedef const E2 expression2_type;
0601
0602
0603 #if BOOST_UBLAS_TYPE_CHECK
0604         typedef typename M::value_type value_type;
0605         matrix<value_type, column_major> cm (m);
0606         typedef typename type_traits<value_type>::real_type real_type;
0607         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
0608         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), column_major_tag ());
0609 #endif
0610         typename expression2_type::const_iterator2 it2 (e2 ().begin2 ());
0611         typename expression2_type::const_iterator2 it2_end (e2 ().end2 ());
0612         while (it2 != it2_end) {
0613 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
0614             typename expression2_type::const_iterator1 it1 (it2.begin ());
0615             typename expression2_type::const_iterator1 it1_end (it2.end ());
0616 #else
0617             typename expression2_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));
0618             typename expression2_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));
0619 #endif
0620             while (it1 != it1_end) {
0621                 // column (m, it2.index2 ()).plus_assign (*it1 * column (e1 (), it1.index1 ()));
0622                 matrix_column<expression1_type> mc (e1 (), it1.index1 ());
0623                 typename matrix_column<expression1_type>::const_iterator itc (mc.begin ());
0624                 typename matrix_column<expression1_type>::const_iterator itc_end (mc.end ());
0625                 while (itc != itc_end) {
0626                     if(triangular_restriction::other (itc.index (), it2.index2 ()))
0627                        m (itc.index (), it2.index2 ()) += *it1 * *itc;
0628                     ++ itc;
0629                 }
0630                 ++ it1;
0631             }
0632             ++ it2;
0633         }
0634 #if BOOST_UBLAS_TYPE_CHECK
0635         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
0636 #endif
0637         return m;
0638     }
0639
0640     // Dispatcher
0641     template<class M, class E1, class E2, class TRI>
0642     BOOST_UBLAS_INLINE
0643     M &
0644     axpy_prod (const matrix_expression<E1> &e1,
0645                const matrix_expression<E2> &e2,
0646                M &m, TRI, bool init = true) {
0647         typedef typename M::value_type value_type;
0648         typedef typename M::storage_category storage_category;
0649         typedef typename M::orientation_category orientation_category;
0650         typedef TRI triangular_restriction;
0651
0652         if (init)
0653             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
0654         return axpy_prod (e1, e2, m, triangular_restriction (), storage_category (), orientation_category ());
0655     }
0656     template<class M, class E1, class E2, class TRI>
0657     BOOST_UBLAS_INLINE
0658     M
0659     axpy_prod (const matrix_expression<E1> &e1,
0660                const matrix_expression<E2> &e2,
0661                TRI) {
0662         typedef M matrix_type;
0663         typedef TRI triangular_restriction;
0664
0665         matrix_type m (e1 ().size1 (), e2 ().size2 ());
0666         return axpy_prod (e1, e2, m, triangular_restriction (), true);
0667     }
0668
0669   /** \brief computes <tt>M += A X</tt> or <tt>M = A X</tt> in an
0670           optimized fashion.
0671
0672           \param e1 the matrix expression \c A
0673           \param e2 the matrix expression \c X
0674           \param m  the result matrix \c M
0675           \param init a boolean parameter
0676
0677           <tt>axpy_prod(A, X, M, init)</tt> implements the well known
0678           axpy-product.  Setting \a init to \c true is equivalent to call
0679           <tt>M.clear()</tt> before <tt>axpy_prod</tt>. Currently \a init
0680           defaults to \c true, but this may change in the future.
0681
0682           Up to now there are no specialisations.
0683
0684           \ingroup blas3
0685
0686           \internal
0687
0688           template parameters:
0689           \param M type of the result matrix \c M
0690           \param E1 type of a matrix expression \c A
0691           \param E2 type of a matrix expression \c X
0692   */
0693     template<class M, class E1, class E2>
0694     BOOST_UBLAS_INLINE
0695     M &
0696     axpy_prod (const matrix_expression<E1> &e1,
0697                const matrix_expression<E2> &e2,
0698                M &m, bool init = true) {
0699         typedef typename M::value_type value_type;
0700         typedef typename M::storage_category storage_category;
0701         typedef typename M::orientation_category orientation_category;
0702
0703         if (init)
0704             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
0705         return axpy_prod (e1, e2, m, full (), storage_category (), orientation_category ());
0706     }
0707     template<class M, class E1, class E2>
0708     BOOST_UBLAS_INLINE
0709     M
0710     axpy_prod (const matrix_expression<E1> &e1,
0711                const matrix_expression<E2> &e2) {
0712         typedef M matrix_type;
0713
0714         matrix_type m (e1 ().size1 (), e2 ().size2 ());
0715         return axpy_prod (e1, e2, m, full (), true);
0716     }
0717
0718
0719     template<class M, class E1, class E2>
0720     BOOST_UBLAS_INLINE
0721     M &
0722     opb_prod (const matrix_expression<E1> &e1,
0723               const matrix_expression<E2> &e2,
0724               M &m,
0725               dense_proxy_tag, row_major_tag) {
0726         typedef typename M::size_type size_type;
0727         typedef typename M::value_type value_type;
0728
0729 #if BOOST_UBLAS_TYPE_CHECK
0730         matrix<value_type, row_major> cm (m);
0731         typedef typename type_traits<value_type>::real_type real_type;
0732         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
0733         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), row_major_tag ());
0734 #endif
0735         size_type size (BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ()));
0736         for (size_type k = 0; k < size; ++ k) {
0737             vector<value_type> ce1 (column (e1 (), k));
0738             vector<value_type> re2 (row (e2 (), k));
0739             m.plus_assign (outer_prod (ce1, re2));
0740         }
0741 #if BOOST_UBLAS_TYPE_CHECK
0742         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
0743 #endif
0744         return m;
0745     }
0746
0747     template<class M, class E1, class E2>
0748     BOOST_UBLAS_INLINE
0749     M &
0750     opb_prod (const matrix_expression<E1> &e1,
0751               const matrix_expression<E2> &e2,
0752               M &m,
0753               dense_proxy_tag, column_major_tag) {
0754         typedef typename M::size_type size_type;
0755         typedef typename M::value_type value_type;
0756
0757 #if BOOST_UBLAS_TYPE_CHECK
0758         matrix<value_type, column_major> cm (m);
0759         typedef typename type_traits<value_type>::real_type real_type;
0760         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
0761         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), column_major_tag ());
0762 #endif
0763         size_type size (BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ()));
0764         for (size_type k = 0; k < size; ++ k) {
0765             vector<value_type> ce1 (column (e1 (), k));
0766             vector<value_type> re2 (row (e2 (), k));
0767             m.plus_assign (outer_prod (ce1, re2));
0768         }
0769 #if BOOST_UBLAS_TYPE_CHECK
0770         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
0771 #endif
0772         return m;
0773     }
0774
0775     // Dispatcher
0776
0777   /** \brief computes <tt>M += A X</tt> or <tt>M = A X</tt> in an
0778           optimized fashion.
0779
0780           \param e1 the matrix expression \c A
0781           \param e2 the matrix expression \c X
0782           \param m  the result matrix \c M
0783           \param init a boolean parameter
0784
0785           <tt>opb_prod(A, X, M, init)</tt> implements the well known
0786           axpy-product. Setting \a init to \c true is equivalent to call
0787           <tt>M.clear()</tt> before <tt>opb_prod</tt>. Currently \a init
0788           defaults to \c true, but this may change in the future.
0789
0790           This function may give a speedup if \c A has less columns than
0791           rows, because the product is computed as a sum of outer
0792           products.
0793
0794           \ingroup blas3
0795
0796           \internal
0797
0798           template parameters:
0799           \param M type of the result matrix \c M
0800           \param E1 type of a matrix expression \c A
0801           \param E2 type of a matrix expression \c X
0802   */
0803     template<class M, class E1, class E2>
0804     BOOST_UBLAS_INLINE
0805     M &
0806     opb_prod (const matrix_expression<E1> &e1,
0807               const matrix_expression<E2> &e2,
0808               M &m, bool init = true) {
0809         typedef typename M::value_type value_type;
0810         typedef typename M::storage_category storage_category;
0811         typedef typename M::orientation_category orientation_category;
0812
0813         if (init)
0814             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
0815         return opb_prod (e1, e2, m, storage_category (), orientation_category ());
0816     }
0817     template<class M, class E1, class E2>
0818     BOOST_UBLAS_INLINE
0819     M
0820     opb_prod (const matrix_expression<E1> &e1,
0821               const matrix_expression<E2> &e2) {
0822         typedef M matrix_type;
0823
0824         matrix_type m (e1 ().size1 (), e2 ().size2 ());
0825         return opb_prod (e1, e2, m, true);
0826     }
0827
0828 }}}
0829
0830 #endif