...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
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