ublas: operation.hpp Source File

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