- Main Page
- Related Pages
- Classes
- Files

- File List

# operation.hpp

    00001 //
    00002 //  Copyright (c) 2000-2002
    00003 //  Joerg Walter, Mathias Koch
    00004 //
    00005 //  Distributed under the Boost Software License, Version 1.0. (See
    00006 //  accompanying file LICENSE_1_0.txt or copy at
    00007 //  http://www.boost.org/LICENSE_1_0.txt)
    00008 //
    00009 //  The authors gratefully acknowledge the support of
    00010 //  GeNeSys mbH & Co. KG in producing this work.
    00011 //
    00012
    00013 #ifndef _BOOST_UBLAS_OPERATION_
    00014 #define _BOOST_UBLAS_OPERATION_
    00015
    00016 #include <boost/numeric/ublas/matrix_proxy.hpp>
    00017
    00022 // axpy-based products
    00023 // Alexei Novakov had a lot of ideas to improve these. Thanks.
    00024 // Hendrik Kueck proposed some new kernel. Thanks again.
    00025
    00026 namespace boost { namespace numeric { namespace ublas {
    00027
    00028     template<class V, class T1, class L1, class IA1, class TA1, class E2>
    00029     BOOST_UBLAS_INLINE
    00030     V &
    00031     axpy_prod (const compressed_matrix<T1, L1, 0, IA1, TA1> &e1,
    00032                const vector_expression<E2> &e2,
    00033                V &v, row_major_tag) {
    00034         typedef typename V::size_type size_type;
    00035         typedef typename V::value_type value_type;
    00036
    00037         for (size_type i = 0; i < e1.filled1 () -1; ++ i) {
    00038             size_type begin = e1.index1_data () [i];
    00039             size_type end = e1.index1_data () [i + 1];
    00040             value_type t (v (i));
    00041             for (size_type j = begin; j < end; ++ j)
    00042                 t += e1.value_data () [j] * e2 () (e1.index2_data () [j]);
    00043             v (i) = t;
    00044         }
    00045         return v;
    00046     }
    00047
    00048     template<class V, class T1, class L1, class IA1, class TA1, class E2>
    00049     BOOST_UBLAS_INLINE
    00050     V &
    00051     axpy_prod (const compressed_matrix<T1, L1, 0, IA1, TA1> &e1,
    00052                const vector_expression<E2> &e2,
    00053                V &v, column_major_tag) {
    00054         typedef typename V::size_type size_type;
    00055
    00056         for (size_type j = 0; j < e1.filled1 () -1; ++ j) {
    00057             size_type begin = e1.index1_data () [j];
    00058             size_type end = e1.index1_data () [j + 1];
    00059             for (size_type i = begin; i < end; ++ i)
    00060                 v (e1.index2_data () [i]) += e1.value_data () [i] * e2 () (j);
    00061         }
    00062         return v;
    00063     }
    00064
    00065     // Dispatcher
    00066     template<class V, class T1, class L1, class IA1, class TA1, class E2>
    00067     BOOST_UBLAS_INLINE
    00068     V &
    00069     axpy_prod (const compressed_matrix<T1, L1, 0, IA1, TA1> &e1,
    00070                const vector_expression<E2> &e2,
    00071                V &v, bool init = true) {
    00072         typedef typename V::value_type value_type;
    00073         typedef typename L1::orientation_category orientation_category;
    00074
    00075         if (init)
    00076             v.assign (zero_vector<value_type> (e1.size1 ()));
    00077 #if BOOST_UBLAS_TYPE_CHECK
    00078         vector<value_type> cv (v);
    00079         typedef typename type_traits<value_type>::real_type real_type;
    00080         real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
    00081         indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));
    00082 #endif
    00083         axpy_prod (e1, e2, v, orientation_category ());
    00084 #if BOOST_UBLAS_TYPE_CHECK
    00085         BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
    00086 #endif
    00087         return v;
    00088     }
    00089     template<class V, class T1, class L1, class IA1, class TA1, class E2>
    00090     BOOST_UBLAS_INLINE
    00091     V
    00092     axpy_prod (const compressed_matrix<T1, L1, 0, IA1, TA1> &e1,
    00093                const vector_expression<E2> &e2) {
    00094         typedef V vector_type;
    00095
    00096         vector_type v (e1.size1 ());
    00097         return axpy_prod (e1, e2, v, true);
    00098     }
    00099
    00100     template<class V, class T1, class L1, class IA1, class TA1, class E2>
    00101     BOOST_UBLAS_INLINE
    00102     V &
    00103     axpy_prod (const coordinate_matrix<T1, L1, 0, IA1, TA1> &e1,
    00104                const vector_expression<E2> &e2,
    00105                V &v, bool init = true) {
    00106         typedef typename V::size_type size_type;
    00107         typedef typename V::value_type value_type;
    00108         typedef L1 layout_type;
    00109
    00110         size_type size1 = e1.size1();
    00111         size_type size2 = e1.size2();
    00112
    00113         if (init) {
    00114             noalias(v) = zero_vector<value_type>(size1);
    00115         }
    00116
    00117         for (size_type i = 0; i < e1.nnz(); ++i) {
    00118             size_type row_index = layout_type::index_M( e1.index1_data () [i], e1.index2_data () [i] );
    00119             size_type col_index = layout_type::index_m( e1.index1_data () [i], e1.index2_data () [i] );
    00120             v( row_index ) += e1.value_data () [i] * e2 () (col_index);
    00121         }
    00122         return v;
    00123     }
    00124
    00125     template<class V, class E1, class E2>
    00126     BOOST_UBLAS_INLINE
    00127     V &
    00128     axpy_prod (const matrix_expression<E1> &e1,
    00129                const vector_expression<E2> &e2,
    00130                V &v, packed_random_access_iterator_tag, row_major_tag) {
    00131         typedef const E1 expression1_type;
    00132         typedef const E2 expression2_type;
    00133         typedef typename V::size_type size_type;
    00134
    00135         typename expression1_type::const_iterator1 it1 (e1 ().begin1 ());
    00136         typename expression1_type::const_iterator1 it1_end (e1 ().end1 ());
    00137         while (it1 != it1_end) {
    00138             size_type index1 (it1.index1 ());
    00139 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    00140             typename expression1_type::const_iterator2 it2 (it1.begin ());
    00141             typename expression1_type::const_iterator2 it2_end (it1.end ());
    00142 #else
    00143             typename expression1_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));
    00144             typename expression1_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));
    00145 #endif
    00146             while (it2 != it2_end) {
    00147                 v (index1) += *it2 * e2 () (it2.index2 ());
    00148                 ++ it2;
    00149             }
    00150             ++ it1;
    00151         }
    00152         return v;
    00153     }
    00154
    00155     template<class V, class E1, class E2>
    00156     BOOST_UBLAS_INLINE
    00157     V &
    00158     axpy_prod (const matrix_expression<E1> &e1,
    00159                const vector_expression<E2> &e2,
    00160                V &v, packed_random_access_iterator_tag, column_major_tag) {
    00161         typedef const E1 expression1_type;
    00162         typedef const E2 expression2_type;
    00163         typedef typename V::size_type size_type;
    00164
    00165         typename expression1_type::const_iterator2 it2 (e1 ().begin2 ());
    00166         typename expression1_type::const_iterator2 it2_end (e1 ().end2 ());
    00167         while (it2 != it2_end) {
    00168             size_type index2 (it2.index2 ());
    00169 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    00170             typename expression1_type::const_iterator1 it1 (it2.begin ());
    00171             typename expression1_type::const_iterator1 it1_end (it2.end ());
    00172 #else
    00173             typename expression1_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));
    00174             typename expression1_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));
    00175 #endif
    00176             while (it1 != it1_end) {
    00177                 v (it1.index1 ()) += *it1 * e2 () (index2);
    00178                 ++ it1;
    00179             }
    00180             ++ it2;
    00181         }
    00182         return v;
    00183     }
    00184
    00185     template<class V, class E1, class E2>
    00186     BOOST_UBLAS_INLINE
    00187     V &
    00188     axpy_prod (const matrix_expression<E1> &e1,
    00189                const vector_expression<E2> &e2,
    00190                V &v, sparse_bidirectional_iterator_tag) {
    00191         typedef const E1 expression1_type;
    00192         typedef const E2 expression2_type;
    00193         typedef typename V::size_type size_type;
    00194
    00195         typename expression2_type::const_iterator it (e2 ().begin ());
    00196         typename expression2_type::const_iterator it_end (e2 ().end ());
    00197         while (it != it_end) {
    00198             v.plus_assign (column (e1 (), it.index ()) * *it);
    00199             ++ it;
    00200         }
    00201         return v;
    00202     }
    00203
    00204     // Dispatcher
    00205     template<class V, class E1, class E2>
    00206     BOOST_UBLAS_INLINE
    00207     V &
    00208     axpy_prod (const matrix_expression<E1> &e1,
    00209                const vector_expression<E2> &e2,
    00210                V &v, packed_random_access_iterator_tag) {
    00211         typedef typename E1::orientation_category orientation_category;
    00212         return axpy_prod (e1, e2, v, packed_random_access_iterator_tag (), orientation_category ());
    00213     }
    00214
    00215
    00241     template<class V, class E1, class E2>
    00242     BOOST_UBLAS_INLINE
    00243     V &
    00244     axpy_prod (const matrix_expression<E1> &e1,
    00245                const vector_expression<E2> &e2,
    00246                V &v, bool init = true) {
    00247         typedef typename V::value_type value_type;
    00248         typedef typename E2::const_iterator::iterator_category iterator_category;
    00249
    00250         if (init)
    00251             v.assign (zero_vector<value_type> (e1 ().size1 ()));
    00252 #if BOOST_UBLAS_TYPE_CHECK
    00253         vector<value_type> cv (v);
    00254         typedef typename type_traits<value_type>::real_type real_type;
    00255         real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
    00256         indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));
    00257 #endif
    00258         axpy_prod (e1, e2, v, iterator_category ());
    00259 #if BOOST_UBLAS_TYPE_CHECK
    00260         BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
    00261 #endif
    00262         return v;
    00263     }
    00264     template<class V, class E1, class E2>
    00265     BOOST_UBLAS_INLINE
    00266     V
    00267     axpy_prod (const matrix_expression<E1> &e1,
    00268                const vector_expression<E2> &e2) {
    00269         typedef V vector_type;
    00270
    00271         vector_type v (e1 ().size1 ());
    00272         return axpy_prod (e1, e2, v, true);
    00273     }
    00274
    00275     template<class V, class E1, class T2, class IA2, class TA2>
    00276     BOOST_UBLAS_INLINE
    00277     V &
    00278     axpy_prod (const vector_expression<E1> &e1,
    00279                const compressed_matrix<T2, column_major, 0, IA2, TA2> &e2,
    00280                V &v, column_major_tag) {
    00281         typedef typename V::size_type size_type;
    00282         typedef typename V::value_type value_type;
    00283
    00284         for (size_type j = 0; j < e2.filled1 () -1; ++ j) {
    00285             size_type begin = e2.index1_data () [j];
    00286             size_type end = e2.index1_data () [j + 1];
    00287             value_type t (v (j));
    00288             for (size_type i = begin; i < end; ++ i)
    00289                 t += e2.value_data () [i] * e1 () (e2.index2_data () [i]);
    00290             v (j) = t;
    00291         }
    00292         return v;
    00293     }
    00294
    00295     template<class V, class E1, class T2, class IA2, class TA2>
    00296     BOOST_UBLAS_INLINE
    00297     V &
    00298     axpy_prod (const vector_expression<E1> &e1,
    00299                const compressed_matrix<T2, row_major, 0, IA2, TA2> &e2,
    00300                V &v, row_major_tag) {
    00301         typedef typename V::size_type size_type;
    00302
    00303         for (size_type i = 0; i < e2.filled1 () -1; ++ i) {
    00304             size_type begin = e2.index1_data () [i];
    00305             size_type end = e2.index1_data () [i + 1];
    00306             for (size_type j = begin; j < end; ++ j)
    00307                 v (e2.index2_data () [j]) += e2.value_data () [j] * e1 () (i);
    00308         }
    00309         return v;
    00310     }
    00311
    00312     // Dispatcher
    00313     template<class V, class E1, class T2, class L2, class IA2, class TA2>
    00314     BOOST_UBLAS_INLINE
    00315     V &
    00316     axpy_prod (const vector_expression<E1> &e1,
    00317                const compressed_matrix<T2, L2, 0, IA2, TA2> &e2,
    00318                V &v, bool init = true) {
    00319         typedef typename V::value_type value_type;
    00320         typedef typename L2::orientation_category orientation_category;
    00321
    00322         if (init)
    00323             v.assign (zero_vector<value_type> (e2.size2 ()));
    00324 #if BOOST_UBLAS_TYPE_CHECK
    00325         vector<value_type> cv (v);
    00326         typedef typename type_traits<value_type>::real_type real_type;
    00327         real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
    00328         indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));
    00329 #endif
    00330         axpy_prod (e1, e2, v, orientation_category ());
    00331 #if BOOST_UBLAS_TYPE_CHECK
    00332         BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
    00333 #endif
    00334         return v;
    00335     }
    00336     template<class V, class E1, class T2, class L2, class IA2, class TA2>
    00337     BOOST_UBLAS_INLINE
    00338     V
    00339     axpy_prod (const vector_expression<E1> &e1,
    00340                const compressed_matrix<T2, L2, 0, IA2, TA2> &e2) {
    00341         typedef V vector_type;
    00342
    00343         vector_type v (e2.size2 ());
    00344         return axpy_prod (e1, e2, v, true);
    00345     }
    00346
    00347     template<class V, class E1, class E2>
    00348     BOOST_UBLAS_INLINE
    00349     V &
    00350     axpy_prod (const vector_expression<E1> &e1,
    00351                const matrix_expression<E2> &e2,
    00352                V &v, packed_random_access_iterator_tag, column_major_tag) {
    00353         typedef const E1 expression1_type;
    00354         typedef const E2 expression2_type;
    00355         typedef typename V::size_type size_type;
    00356
    00357         typename expression2_type::const_iterator2 it2 (e2 ().begin2 ());
    00358         typename expression2_type::const_iterator2 it2_end (e2 ().end2 ());
    00359         while (it2 != it2_end) {
    00360             size_type index2 (it2.index2 ());
    00361 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    00362             typename expression2_type::const_iterator1 it1 (it2.begin ());
    00363             typename expression2_type::const_iterator1 it1_end (it2.end ());
    00364 #else
    00365             typename expression2_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));
    00366             typename expression2_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));
    00367 #endif
    00368             while (it1 != it1_end) {
    00369                 v (index2) += *it1 * e1 () (it1.index1 ());
    00370                 ++ it1;
    00371             }
    00372             ++ it2;
    00373         }
    00374         return v;
    00375     }
    00376
    00377     template<class V, class E1, class E2>
    00378     BOOST_UBLAS_INLINE
    00379     V &
    00380     axpy_prod (const vector_expression<E1> &e1,
    00381                const matrix_expression<E2> &e2,
    00382                V &v, packed_random_access_iterator_tag, row_major_tag) {
    00383         typedef const E1 expression1_type;
    00384         typedef const E2 expression2_type;
    00385         typedef typename V::size_type size_type;
    00386
    00387         typename expression2_type::const_iterator1 it1 (e2 ().begin1 ());
    00388         typename expression2_type::const_iterator1 it1_end (e2 ().end1 ());
    00389         while (it1 != it1_end) {
    00390             size_type index1 (it1.index1 ());
    00391 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    00392             typename expression2_type::const_iterator2 it2 (it1.begin ());
    00393             typename expression2_type::const_iterator2 it2_end (it1.end ());
    00394 #else
    00395             typename expression2_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));
    00396             typename expression2_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));
    00397 #endif
    00398             while (it2 != it2_end) {
    00399                 v (it2.index2 ()) += *it2 * e1 () (index1);
    00400                 ++ it2;
    00401             }
    00402             ++ it1;
    00403         }
    00404         return v;
    00405     }
    00406
    00407     template<class V, class E1, class E2>
    00408     BOOST_UBLAS_INLINE
    00409     V &
    00410     axpy_prod (const vector_expression<E1> &e1,
    00411                const matrix_expression<E2> &e2,
    00412                V &v, sparse_bidirectional_iterator_tag) {
    00413         typedef const E1 expression1_type;
    00414         typedef const E2 expression2_type;
    00415         typedef typename V::size_type size_type;
    00416
    00417         typename expression1_type::const_iterator it (e1 ().begin ());
    00418         typename expression1_type::const_iterator it_end (e1 ().end ());
    00419         while (it != it_end) {
    00420             v.plus_assign (*it * row (e2 (), it.index ()));
    00421             ++ it;
    00422         }
    00423         return v;
    00424     }
    00425
    00426     // Dispatcher
    00427     template<class V, class E1, class E2>
    00428     BOOST_UBLAS_INLINE
    00429     V &
    00430     axpy_prod (const vector_expression<E1> &e1,
    00431                const matrix_expression<E2> &e2,
    00432                V &v, packed_random_access_iterator_tag) {
    00433         typedef typename E2::orientation_category orientation_category;
    00434         return axpy_prod (e1, e2, v, packed_random_access_iterator_tag (), orientation_category ());
    00435     }
    00436
    00437
    00463     template<class V, class E1, class E2>
    00464     BOOST_UBLAS_INLINE
    00465     V &
    00466     axpy_prod (const vector_expression<E1> &e1,
    00467                const matrix_expression<E2> &e2,
    00468                V &v, bool init = true) {
    00469         typedef typename V::value_type value_type;
    00470         typedef typename E1::const_iterator::iterator_category iterator_category;
    00471
    00472         if (init)
    00473             v.assign (zero_vector<value_type> (e2 ().size2 ()));
    00474 #if BOOST_UBLAS_TYPE_CHECK
    00475         vector<value_type> cv (v);
    00476         typedef typename type_traits<value_type>::real_type real_type;
    00477         real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
    00478         indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));
    00479 #endif
    00480         axpy_prod (e1, e2, v, iterator_category ());
    00481 #if BOOST_UBLAS_TYPE_CHECK
    00482         BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
    00483 #endif
    00484         return v;
    00485     }
    00486     template<class V, class E1, class E2>
    00487     BOOST_UBLAS_INLINE
    00488     V
    00489     axpy_prod (const vector_expression<E1> &e1,
    00490                const matrix_expression<E2> &e2) {
    00491         typedef V vector_type;
    00492
    00493         vector_type v (e2 ().size2 ());
    00494         return axpy_prod (e1, e2, v, true);
    00495     }
    00496
    00497     template<class M, class E1, class E2, class TRI>
    00498     BOOST_UBLAS_INLINE
    00499     M &
    00500     axpy_prod (const matrix_expression<E1> &e1,
    00501                const matrix_expression<E2> &e2,
    00502                M &m, TRI,
    00503                dense_proxy_tag, row_major_tag) {
    00504         typedef M matrix_type;
    00505         typedef const E1 expression1_type;
    00506         typedef const E2 expression2_type;
    00507         typedef typename M::size_type size_type;
    00508         typedef typename M::value_type value_type;
    00509
    00510 #if BOOST_UBLAS_TYPE_CHECK
    00511         matrix<value_type, row_major> cm (m);
    00512         typedef typename type_traits<value_type>::real_type real_type;
    00513         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
    00514         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), row_major_tag ());
    00515 #endif
    00516         size_type size1 (e1 ().size1 ());
    00517         size_type size2 (e1 ().size2 ());
    00518         for (size_type i = 0; i < size1; ++ i)
    00519             for (size_type j = 0; j < size2; ++ j)
    00520                 row (m, i).plus_assign (e1 () (i, j) * row (e2 (), j));
    00521 #if BOOST_UBLAS_TYPE_CHECK
    00522         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
    00523 #endif
    00524         return m;
    00525     }
    00526     template<class M, class E1, class E2, class TRI>
    00527     BOOST_UBLAS_INLINE
    00528     M &
    00529     axpy_prod (const matrix_expression<E1> &e1,
    00530                const matrix_expression<E2> &e2,
    00531                M &m, TRI,
    00532                sparse_proxy_tag, row_major_tag) {
    00533         typedef M matrix_type;
    00534         typedef TRI triangular_restriction;
    00535         typedef const E1 expression1_type;
    00536         typedef const E2 expression2_type;
    00537         typedef typename M::size_type size_type;
    00538         typedef typename M::value_type value_type;
    00539
    00540 #if BOOST_UBLAS_TYPE_CHECK
    00541         matrix<value_type, row_major> cm (m);
    00542         typedef typename type_traits<value_type>::real_type real_type;
    00543         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
    00544         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), row_major_tag ());
    00545 #endif
    00546         typename expression1_type::const_iterator1 it1 (e1 ().begin1 ());
    00547         typename expression1_type::const_iterator1 it1_end (e1 ().end1 ());
    00548         while (it1 != it1_end) {
    00549 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    00550             typename expression1_type::const_iterator2 it2 (it1.begin ());
    00551             typename expression1_type::const_iterator2 it2_end (it1.end ());
    00552 #else
    00553             typename expression1_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));
    00554             typename expression1_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));
    00555 #endif
    00556             while (it2 != it2_end) {
    00557                 // row (m, it1.index1 ()).plus_assign (*it2 * row (e2 (), it2.index2 ()));
    00558                 matrix_row<expression2_type> mr (e2 (), it2.index2 ());
    00559                 typename matrix_row<expression2_type>::const_iterator itr (mr.begin ());
    00560                 typename matrix_row<expression2_type>::const_iterator itr_end (mr.end ());
    00561                 while (itr != itr_end) {
    00562                     if (triangular_restriction::other (it1.index1 (), itr.index ()))
    00563                         m (it1.index1 (), itr.index ()) += *it2 * *itr;
    00564                     ++ itr;
    00565                 }
    00566                 ++ it2;
    00567             }
    00568             ++ it1;
    00569         }
    00570 #if BOOST_UBLAS_TYPE_CHECK
    00571         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
    00572 #endif
    00573         return m;
    00574     }
    00575
    00576     template<class M, class E1, class E2, class TRI>
    00577     BOOST_UBLAS_INLINE
    00578     M &
    00579     axpy_prod (const matrix_expression<E1> &e1,
    00580                const matrix_expression<E2> &e2,
    00581                M &m, TRI,
    00582                dense_proxy_tag, column_major_tag) {
    00583         typedef M matrix_type;
    00584         typedef const E1 expression1_type;
    00585         typedef const E2 expression2_type;
    00586         typedef typename M::size_type size_type;
    00587         typedef typename M::value_type value_type;
    00588
    00589 #if BOOST_UBLAS_TYPE_CHECK
    00590         matrix<value_type, column_major> cm (m);
    00591         typedef typename type_traits<value_type>::real_type real_type;
    00592         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
    00593         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), column_major_tag ());
    00594 #endif
    00595         size_type size1 (e2 ().size1 ());
    00596         size_type size2 (e2 ().size2 ());
    00597         for (size_type j = 0; j < size2; ++ j)
    00598             for (size_type i = 0; i < size1; ++ i)
    00599                 column (m, j).plus_assign (e2 () (i, j) * column (e1 (), i));
    00600 #if BOOST_UBLAS_TYPE_CHECK
    00601         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
    00602 #endif
    00603         return m;
    00604     }
    00605     template<class M, class E1, class E2, class TRI>
    00606     BOOST_UBLAS_INLINE
    00607     M &
    00608     axpy_prod (const matrix_expression<E1> &e1,
    00609                const matrix_expression<E2> &e2,
    00610                M &m, TRI,
    00611                sparse_proxy_tag, column_major_tag) {
    00612         typedef M matrix_type;
    00613         typedef TRI triangular_restriction;
    00614         typedef const E1 expression1_type;
    00615         typedef const E2 expression2_type;
    00616         typedef typename M::size_type size_type;
    00617         typedef typename M::value_type value_type;
    00618
    00619 #if BOOST_UBLAS_TYPE_CHECK
    00620         matrix<value_type, column_major> cm (m);
    00621         typedef typename type_traits<value_type>::real_type real_type;
    00622         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
    00623         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), column_major_tag ());
    00624 #endif
    00625         typename expression2_type::const_iterator2 it2 (e2 ().begin2 ());
    00626         typename expression2_type::const_iterator2 it2_end (e2 ().end2 ());
    00627         while (it2 != it2_end) {
    00628 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    00629             typename expression2_type::const_iterator1 it1 (it2.begin ());
    00630             typename expression2_type::const_iterator1 it1_end (it2.end ());
    00631 #else
    00632             typename expression2_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));
    00633             typename expression2_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));
    00634 #endif
    00635             while (it1 != it1_end) {
    00636                 // column (m, it2.index2 ()).plus_assign (*it1 * column (e1 (), it1.index1 ()));
    00637                 matrix_column<expression1_type> mc (e1 (), it1.index1 ());
    00638                 typename matrix_column<expression1_type>::const_iterator itc (mc.begin ());
    00639                 typename matrix_column<expression1_type>::const_iterator itc_end (mc.end ());
    00640                 while (itc != itc_end) {
    00641                     if(triangular_restriction::other (itc.index (), it2.index2 ()))
    00642                        m (itc.index (), it2.index2 ()) += *it1 * *itc;
    00643                     ++ itc;
    00644                 }
    00645                 ++ it1;
    00646             }
    00647             ++ it2;
    00648         }
    00649 #if BOOST_UBLAS_TYPE_CHECK
    00650         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
    00651 #endif
    00652         return m;
    00653     }
    00654
    00655     // Dispatcher
    00656     template<class M, class E1, class E2, class TRI>
    00657     BOOST_UBLAS_INLINE
    00658     M &
    00659     axpy_prod (const matrix_expression<E1> &e1,
    00660                const matrix_expression<E2> &e2,
    00661                M &m, TRI, bool init = true) {
    00662         typedef typename M::value_type value_type;
    00663         typedef typename M::storage_category storage_category;
    00664         typedef typename M::orientation_category orientation_category;
    00665         typedef TRI triangular_restriction;
    00666
    00667         if (init)
    00668             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
    00669         return axpy_prod (e1, e2, m, triangular_restriction (), storage_category (), orientation_category ());
    00670     }
    00671     template<class M, class E1, class E2, class TRI>
    00672     BOOST_UBLAS_INLINE
    00673     M
    00674     axpy_prod (const matrix_expression<E1> &e1,
    00675                const matrix_expression<E2> &e2,
    00676                TRI) {
    00677         typedef M matrix_type;
    00678         typedef TRI triangular_restriction;
    00679
    00680         matrix_type m (e1 ().size1 (), e2 ().size2 ());
    00681         return axpy_prod (e1, e2, m, triangular_restriction (), true);
    00682     }
    00683
    00708     template<class M, class E1, class E2>
    00709     BOOST_UBLAS_INLINE
    00710     M &
    00711     axpy_prod (const matrix_expression<E1> &e1,
    00712                const matrix_expression<E2> &e2,
    00713                M &m, bool init = true) {
    00714         typedef typename M::value_type value_type;
    00715         typedef typename M::storage_category storage_category;
    00716         typedef typename M::orientation_category orientation_category;
    00717
    00718         if (init)
    00719             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
    00720         return axpy_prod (e1, e2, m, full (), storage_category (), orientation_category ());
    00721     }
    00722     template<class M, class E1, class E2>
    00723     BOOST_UBLAS_INLINE
    00724     M
    00725     axpy_prod (const matrix_expression<E1> &e1,
    00726                const matrix_expression<E2> &e2) {
    00727         typedef M matrix_type;
    00728
    00729         matrix_type m (e1 ().size1 (), e2 ().size2 ());
    00730         return axpy_prod (e1, e2, m, full (), true);
    00731     }
    00732
    00733
    00734     template<class M, class E1, class E2>
    00735     BOOST_UBLAS_INLINE
    00736     M &
    00737     opb_prod (const matrix_expression<E1> &e1,
    00738               const matrix_expression<E2> &e2,
    00739               M &m,
    00740               dense_proxy_tag, row_major_tag) {
    00741         typedef M matrix_type;
    00742         typedef const E1 expression1_type;
    00743         typedef const E2 expression2_type;
    00744         typedef typename M::size_type size_type;
    00745         typedef typename M::value_type value_type;
    00746
    00747 #if BOOST_UBLAS_TYPE_CHECK
    00748         matrix<value_type, row_major> cm (m);
    00749         typedef typename type_traits<value_type>::real_type real_type;
    00750         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
    00751         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), row_major_tag ());
    00752 #endif
    00753         size_type size (BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ()));
    00754         for (size_type k = 0; k < size; ++ k) {
    00755             vector<value_type> ce1 (column (e1 (), k));
    00756             vector<value_type> re2 (row (e2 (), k));
    00757             m.plus_assign (outer_prod (ce1, re2));
    00758         }
    00759 #if BOOST_UBLAS_TYPE_CHECK
    00760         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
    00761 #endif
    00762         return m;
    00763     }
    00764
    00765     template<class M, class E1, class E2>
    00766     BOOST_UBLAS_INLINE
    00767     M &
    00768     opb_prod (const matrix_expression<E1> &e1,
    00769               const matrix_expression<E2> &e2,
    00770               M &m,
    00771               dense_proxy_tag, column_major_tag) {
    00772         typedef M matrix_type;
    00773         typedef const E1 expression1_type;
    00774         typedef const E2 expression2_type;
    00775         typedef typename M::size_type size_type;
    00776         typedef typename M::value_type value_type;
    00777
    00778 #if BOOST_UBLAS_TYPE_CHECK
    00779         matrix<value_type, column_major> cm (m);
    00780         typedef typename type_traits<value_type>::real_type real_type;
    00781         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
    00782         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), column_major_tag ());
    00783 #endif
    00784         size_type size (BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ()));
    00785         for (size_type k = 0; k < size; ++ k) {
    00786             vector<value_type> ce1 (column (e1 (), k));
    00787             vector<value_type> re2 (row (e2 (), k));
    00788             m.plus_assign (outer_prod (ce1, re2));
    00789         }
    00790 #if BOOST_UBLAS_TYPE_CHECK
    00791         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
    00792 #endif
    00793         return m;
    00794     }
    00795
    00796     // Dispatcher
    00797
    00824     template<class M, class E1, class E2>
    00825     BOOST_UBLAS_INLINE
    00826     M &
    00827     opb_prod (const matrix_expression<E1> &e1,
    00828               const matrix_expression<E2> &e2,
    00829               M &m, bool init = true) {
    00830         typedef typename M::value_type value_type;
    00831         typedef typename M::storage_category storage_category;
    00832         typedef typename M::orientation_category orientation_category;
    00833
    00834         if (init)
    00835             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
    00836         return opb_prod (e1, e2, m, storage_category (), orientation_category ());
    00837     }
    00838     template<class M, class E1, class E2>
    00839     BOOST_UBLAS_INLINE
    00840     M
    00841     opb_prod (const matrix_expression<E1> &e1,
    00842               const matrix_expression<E2> &e2) {
    00843         typedef M matrix_type;
    00844
    00845         matrix_type m (e1 ().size1 (), e2 ().size2 ());
    00846         return opb_prod (e1, e2, m, true);
    00847     }
    00848
    00849 }}}
    00850
    00851 #endif