...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-2009 00003 // Joerg Walter, Mathias Koch, Gunter Winkler 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_FUNCTIONAL_ 00014 #define _BOOST_UBLAS_FUNCTIONAL_ 00015 00016 #include <functional> 00017 00018 #include <boost/numeric/ublas/traits.hpp> 00019 #ifdef BOOST_UBLAS_USE_DUFF_DEVICE 00020 #include <boost/numeric/ublas/detail/duff.hpp> 00021 #endif 00022 #ifdef BOOST_UBLAS_USE_SIMD 00023 #include <boost/numeric/ublas/detail/raw.hpp> 00024 #else 00025 namespace boost { namespace numeric { namespace ublas { namespace raw { 00026 }}}} 00027 #endif 00028 #ifdef BOOST_UBLAS_HAVE_BINDINGS 00029 #include <boost/numeric/bindings/traits/std_vector.hpp> 00030 #include <boost/numeric/bindings/traits/ublas_vector.hpp> 00031 #include <boost/numeric/bindings/traits/ublas_matrix.hpp> 00032 #include <boost/numeric/bindings/atlas/cblas.hpp> 00033 #endif 00034 00035 #include <boost/numeric/ublas/detail/definitions.hpp> 00036 00037 00038 00039 namespace boost { namespace numeric { namespace ublas { 00040 00041 // Scalar functors 00042 00043 // Unary 00044 template<class T> 00045 struct scalar_unary_functor { 00046 typedef T value_type; 00047 typedef typename type_traits<T>::const_reference argument_type; 00048 typedef typename type_traits<T>::value_type result_type; 00049 }; 00050 00051 template<class T> 00052 struct scalar_identity: 00053 public scalar_unary_functor<T> { 00054 typedef typename scalar_unary_functor<T>::argument_type argument_type; 00055 typedef typename scalar_unary_functor<T>::result_type result_type; 00056 00057 static BOOST_UBLAS_INLINE 00058 result_type apply (argument_type t) { 00059 return t; 00060 } 00061 }; 00062 template<class T> 00063 struct scalar_negate: 00064 public scalar_unary_functor<T> { 00065 typedef typename scalar_unary_functor<T>::argument_type argument_type; 00066 typedef typename scalar_unary_functor<T>::result_type result_type; 00067 00068 static BOOST_UBLAS_INLINE 00069 result_type apply (argument_type t) { 00070 return - t; 00071 } 00072 }; 00073 template<class T> 00074 struct scalar_conj: 00075 public scalar_unary_functor<T> { 00076 typedef typename scalar_unary_functor<T>::value_type value_type; 00077 typedef typename scalar_unary_functor<T>::argument_type argument_type; 00078 typedef typename scalar_unary_functor<T>::result_type result_type; 00079 00080 static BOOST_UBLAS_INLINE 00081 result_type apply (argument_type t) { 00082 return type_traits<value_type>::conj (t); 00083 } 00084 }; 00085 00086 // Unary returning real 00087 template<class T> 00088 struct scalar_real_unary_functor { 00089 typedef T value_type; 00090 typedef typename type_traits<T>::const_reference argument_type; 00091 typedef typename type_traits<T>::real_type result_type; 00092 }; 00093 00094 template<class T> 00095 struct scalar_real: 00096 public scalar_real_unary_functor<T> { 00097 typedef typename scalar_real_unary_functor<T>::value_type value_type; 00098 typedef typename scalar_real_unary_functor<T>::argument_type argument_type; 00099 typedef typename scalar_real_unary_functor<T>::result_type result_type; 00100 00101 static BOOST_UBLAS_INLINE 00102 result_type apply (argument_type t) { 00103 return type_traits<value_type>::real (t); 00104 } 00105 }; 00106 template<class T> 00107 struct scalar_imag: 00108 public scalar_real_unary_functor<T> { 00109 typedef typename scalar_real_unary_functor<T>::value_type value_type; 00110 typedef typename scalar_real_unary_functor<T>::argument_type argument_type; 00111 typedef typename scalar_real_unary_functor<T>::result_type result_type; 00112 00113 static BOOST_UBLAS_INLINE 00114 result_type apply (argument_type t) { 00115 return type_traits<value_type>::imag (t); 00116 } 00117 }; 00118 00119 // Binary 00120 template<class T1, class T2> 00121 struct scalar_binary_functor { 00122 typedef typename type_traits<T1>::const_reference argument1_type; 00123 typedef typename type_traits<T2>::const_reference argument2_type; 00124 typedef typename promote_traits<T1, T2>::promote_type result_type; 00125 }; 00126 00127 template<class T1, class T2> 00128 struct scalar_plus: 00129 public scalar_binary_functor<T1, T2> { 00130 typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type; 00131 typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type; 00132 typedef typename scalar_binary_functor<T1, T2>::result_type result_type; 00133 00134 static BOOST_UBLAS_INLINE 00135 result_type apply (argument1_type t1, argument2_type t2) { 00136 return t1 + t2; 00137 } 00138 }; 00139 template<class T1, class T2> 00140 struct scalar_minus: 00141 public scalar_binary_functor<T1, T2> { 00142 typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type; 00143 typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type; 00144 typedef typename scalar_binary_functor<T1, T2>::result_type result_type; 00145 00146 static BOOST_UBLAS_INLINE 00147 result_type apply (argument1_type t1, argument2_type t2) { 00148 return t1 - t2; 00149 } 00150 }; 00151 template<class T1, class T2> 00152 struct scalar_multiplies: 00153 public scalar_binary_functor<T1, T2> { 00154 typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type; 00155 typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type; 00156 typedef typename scalar_binary_functor<T1, T2>::result_type result_type; 00157 00158 static BOOST_UBLAS_INLINE 00159 result_type apply (argument1_type t1, argument2_type t2) { 00160 return t1 * t2; 00161 } 00162 }; 00163 template<class T1, class T2> 00164 struct scalar_divides: 00165 public scalar_binary_functor<T1, T2> { 00166 typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type; 00167 typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type; 00168 typedef typename scalar_binary_functor<T1, T2>::result_type result_type; 00169 00170 static BOOST_UBLAS_INLINE 00171 result_type apply (argument1_type t1, argument2_type t2) { 00172 return t1 / t2; 00173 } 00174 }; 00175 00176 template<class T1, class T2> 00177 struct scalar_binary_assign_functor { 00178 // ISSUE Remove reference to avoid reference to reference problems 00179 typedef typename type_traits<typename boost::remove_reference<T1>::type>::reference argument1_type; 00180 typedef typename type_traits<T2>::const_reference argument2_type; 00181 }; 00182 00183 struct assign_tag {}; 00184 struct computed_assign_tag {}; 00185 00186 template<class T1, class T2> 00187 struct scalar_assign: 00188 public scalar_binary_assign_functor<T1, T2> { 00189 typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type; 00190 typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type; 00191 #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) 00192 static const bool computed ; 00193 #else 00194 static const bool computed = false ; 00195 #endif 00196 00197 static BOOST_UBLAS_INLINE 00198 void apply (argument1_type t1, argument2_type t2) { 00199 t1 = t2; 00200 } 00201 00202 template<class U1, class U2> 00203 struct rebind { 00204 typedef scalar_assign<U1, U2> other; 00205 }; 00206 }; 00207 00208 #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) 00209 template<class T1, class T2> 00210 const bool scalar_assign<T1,T2>::computed = false; 00211 #endif 00212 00213 template<class T1, class T2> 00214 struct scalar_plus_assign: 00215 public scalar_binary_assign_functor<T1, T2> { 00216 typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type; 00217 typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type; 00218 #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) 00219 static const bool computed ; 00220 #else 00221 static const bool computed = true ; 00222 #endif 00223 00224 static BOOST_UBLAS_INLINE 00225 void apply (argument1_type t1, argument2_type t2) { 00226 t1 += t2; 00227 } 00228 00229 template<class U1, class U2> 00230 struct rebind { 00231 typedef scalar_plus_assign<U1, U2> other; 00232 }; 00233 }; 00234 00235 #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) 00236 template<class T1, class T2> 00237 const bool scalar_plus_assign<T1,T2>::computed = true; 00238 #endif 00239 00240 template<class T1, class T2> 00241 struct scalar_minus_assign: 00242 public scalar_binary_assign_functor<T1, T2> { 00243 typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type; 00244 typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type; 00245 #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) 00246 static const bool computed ; 00247 #else 00248 static const bool computed = true ; 00249 #endif 00250 00251 static BOOST_UBLAS_INLINE 00252 void apply (argument1_type t1, argument2_type t2) { 00253 t1 -= t2; 00254 } 00255 00256 template<class U1, class U2> 00257 struct rebind { 00258 typedef scalar_minus_assign<U1, U2> other; 00259 }; 00260 }; 00261 00262 #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) 00263 template<class T1, class T2> 00264 const bool scalar_minus_assign<T1,T2>::computed = true; 00265 #endif 00266 00267 template<class T1, class T2> 00268 struct scalar_multiplies_assign: 00269 public scalar_binary_assign_functor<T1, T2> { 00270 typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type; 00271 typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type; 00272 static const bool computed = true; 00273 00274 static BOOST_UBLAS_INLINE 00275 void apply (argument1_type t1, argument2_type t2) { 00276 t1 *= t2; 00277 } 00278 00279 template<class U1, class U2> 00280 struct rebind { 00281 typedef scalar_multiplies_assign<U1, U2> other; 00282 }; 00283 }; 00284 template<class T1, class T2> 00285 struct scalar_divides_assign: 00286 public scalar_binary_assign_functor<T1, T2> { 00287 typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type; 00288 typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type; 00289 static const bool computed ; 00290 00291 static BOOST_UBLAS_INLINE 00292 void apply (argument1_type t1, argument2_type t2) { 00293 t1 /= t2; 00294 } 00295 00296 template<class U1, class U2> 00297 struct rebind { 00298 typedef scalar_divides_assign<U1, U2> other; 00299 }; 00300 }; 00301 template<class T1, class T2> 00302 const bool scalar_divides_assign<T1,T2>::computed = true; 00303 00304 template<class T1, class T2> 00305 struct scalar_binary_swap_functor { 00306 typedef typename type_traits<typename boost::remove_reference<T1>::type>::reference argument1_type; 00307 typedef typename type_traits<typename boost::remove_reference<T2>::type>::reference argument2_type; 00308 }; 00309 00310 template<class T1, class T2> 00311 struct scalar_swap: 00312 public scalar_binary_swap_functor<T1, T2> { 00313 typedef typename scalar_binary_swap_functor<T1, T2>::argument1_type argument1_type; 00314 typedef typename scalar_binary_swap_functor<T1, T2>::argument2_type argument2_type; 00315 00316 static BOOST_UBLAS_INLINE 00317 void apply (argument1_type t1, argument2_type t2) { 00318 std::swap (t1, t2); 00319 } 00320 00321 template<class U1, class U2> 00322 struct rebind { 00323 typedef scalar_swap<U1, U2> other; 00324 }; 00325 }; 00326 00327 // Vector functors 00328 00329 // Unary returning scalar 00330 template<class V> 00331 struct vector_scalar_unary_functor { 00332 typedef typename V::value_type value_type; 00333 typedef typename V::value_type result_type; 00334 }; 00335 00336 template<class V> 00337 struct vector_sum: 00338 public vector_scalar_unary_functor<V> { 00339 typedef typename vector_scalar_unary_functor<V>::value_type value_type; 00340 typedef typename vector_scalar_unary_functor<V>::result_type result_type; 00341 00342 template<class E> 00343 static BOOST_UBLAS_INLINE 00344 result_type apply (const vector_expression<E> &e) { 00345 result_type t = result_type (0); 00346 typedef typename E::size_type vector_size_type; 00347 vector_size_type size (e ().size ()); 00348 for (vector_size_type i = 0; i < size; ++ i) 00349 t += e () (i); 00350 return t; 00351 } 00352 // Dense case 00353 template<class D, class I> 00354 static BOOST_UBLAS_INLINE 00355 result_type apply (D size, I it) { 00356 result_type t = result_type (0); 00357 while (-- size >= 0) 00358 t += *it, ++ it; 00359 return t; 00360 } 00361 // Sparse case 00362 template<class I> 00363 static BOOST_UBLAS_INLINE 00364 result_type apply (I it, const I &it_end) { 00365 result_type t = result_type (0); 00366 while (it != it_end) 00367 t += *it, ++ it; 00368 return t; 00369 } 00370 }; 00371 00372 // Unary returning real scalar 00373 template<class V> 00374 struct vector_scalar_real_unary_functor { 00375 typedef typename V::value_type value_type; 00376 typedef typename type_traits<value_type>::real_type real_type; 00377 typedef real_type result_type; 00378 }; 00379 00380 template<class V> 00381 struct vector_norm_1: 00382 public vector_scalar_real_unary_functor<V> { 00383 typedef typename vector_scalar_real_unary_functor<V>::value_type value_type; 00384 typedef typename vector_scalar_real_unary_functor<V>::real_type real_type; 00385 typedef typename vector_scalar_real_unary_functor<V>::result_type result_type; 00386 00387 template<class E> 00388 static BOOST_UBLAS_INLINE 00389 result_type apply (const vector_expression<E> &e) { 00390 real_type t = real_type (); 00391 typedef typename E::size_type vector_size_type; 00392 vector_size_type size (e ().size ()); 00393 for (vector_size_type i = 0; i < size; ++ i) { 00394 real_type u (type_traits<value_type>::type_abs (e () (i))); 00395 t += u; 00396 } 00397 return t; 00398 } 00399 // Dense case 00400 template<class D, class I> 00401 static BOOST_UBLAS_INLINE 00402 result_type apply (D size, I it) { 00403 real_type t = real_type (); 00404 while (-- size >= 0) { 00405 real_type u (type_traits<value_type>::norm_1 (*it)); 00406 t += u; 00407 ++ it; 00408 } 00409 return t; 00410 } 00411 // Sparse case 00412 template<class I> 00413 static BOOST_UBLAS_INLINE 00414 result_type apply (I it, const I &it_end) { 00415 real_type t = real_type (); 00416 while (it != it_end) { 00417 real_type u (type_traits<value_type>::norm_1 (*it)); 00418 t += u; 00419 ++ it; 00420 } 00421 return t; 00422 } 00423 }; 00424 template<class V> 00425 struct vector_norm_2: 00426 public vector_scalar_real_unary_functor<V> { 00427 typedef typename vector_scalar_real_unary_functor<V>::value_type value_type; 00428 typedef typename vector_scalar_real_unary_functor<V>::real_type real_type; 00429 typedef typename vector_scalar_real_unary_functor<V>::result_type result_type; 00430 00431 template<class E> 00432 static BOOST_UBLAS_INLINE 00433 result_type apply (const vector_expression<E> &e) { 00434 #ifndef BOOST_UBLAS_SCALED_NORM 00435 real_type t = real_type (); 00436 typedef typename E::size_type vector_size_type; 00437 vector_size_type size (e ().size ()); 00438 for (vector_size_type i = 0; i < size; ++ i) { 00439 real_type u (type_traits<value_type>::norm_2 (e () (i))); 00440 t += u * u; 00441 } 00442 return type_traits<real_type>::type_sqrt (t); 00443 #else 00444 real_type scale = real_type (); 00445 real_type sum_squares (1); 00446 size_type size (e ().size ()); 00447 for (size_type i = 0; i < size; ++ i) { 00448 real_type u (type_traits<value_type>::norm_2 (e () (i))); 00449 if ( real_type () /* zero */ == u ) continue; 00450 if (scale < u) { 00451 real_type v (scale / u); 00452 sum_squares = sum_squares * v * v + real_type (1); 00453 scale = u; 00454 } else { 00455 real_type v (u / scale); 00456 sum_squares += v * v; 00457 } 00458 } 00459 return scale * type_traits<real_type>::type_sqrt (sum_squares); 00460 #endif 00461 } 00462 // Dense case 00463 template<class D, class I> 00464 static BOOST_UBLAS_INLINE 00465 result_type apply (D size, I it) { 00466 #ifndef BOOST_UBLAS_SCALED_NORM 00467 real_type t = real_type (); 00468 while (-- size >= 0) { 00469 real_type u (type_traits<value_type>::norm_2 (*it)); 00470 t += u * u; 00471 ++ it; 00472 } 00473 return type_traits<real_type>::type_sqrt (t); 00474 #else 00475 real_type scale = real_type (); 00476 real_type sum_squares (1); 00477 while (-- size >= 0) { 00478 real_type u (type_traits<value_type>::norm_2 (*it)); 00479 if (scale < u) { 00480 real_type v (scale / u); 00481 sum_squares = sum_squares * v * v + real_type (1); 00482 scale = u; 00483 } else { 00484 real_type v (u / scale); 00485 sum_squares += v * v; 00486 } 00487 ++ it; 00488 } 00489 return scale * type_traits<real_type>::type_sqrt (sum_squares); 00490 #endif 00491 } 00492 // Sparse case 00493 template<class I> 00494 static BOOST_UBLAS_INLINE 00495 result_type apply (I it, const I &it_end) { 00496 #ifndef BOOST_UBLAS_SCALED_NORM 00497 real_type t = real_type (); 00498 while (it != it_end) { 00499 real_type u (type_traits<value_type>::norm_2 (*it)); 00500 t += u * u; 00501 ++ it; 00502 } 00503 return type_traits<real_type>::type_sqrt (t); 00504 #else 00505 real_type scale = real_type (); 00506 real_type sum_squares (1); 00507 while (it != it_end) { 00508 real_type u (type_traits<value_type>::norm_2 (*it)); 00509 if (scale < u) { 00510 real_type v (scale / u); 00511 sum_squares = sum_squares * v * v + real_type (1); 00512 scale = u; 00513 } else { 00514 real_type v (u / scale); 00515 sum_squares += v * v; 00516 } 00517 ++ it; 00518 } 00519 return scale * type_traits<real_type>::type_sqrt (sum_squares); 00520 #endif 00521 } 00522 }; 00523 template<class V> 00524 struct vector_norm_inf: 00525 public vector_scalar_real_unary_functor<V> { 00526 typedef typename vector_scalar_real_unary_functor<V>::value_type value_type; 00527 typedef typename vector_scalar_real_unary_functor<V>::real_type real_type; 00528 typedef typename vector_scalar_real_unary_functor<V>::result_type result_type; 00529 00530 template<class E> 00531 static BOOST_UBLAS_INLINE 00532 result_type apply (const vector_expression<E> &e) { 00533 real_type t = real_type (); 00534 typedef typename E::size_type vector_size_type; 00535 vector_size_type size (e ().size ()); 00536 for (vector_size_type i = 0; i < size; ++ i) { 00537 real_type u (type_traits<value_type>::norm_inf (e () (i))); 00538 if (u > t) 00539 t = u; 00540 } 00541 return t; 00542 } 00543 // Dense case 00544 template<class D, class I> 00545 static BOOST_UBLAS_INLINE 00546 result_type apply (D size, I it) { 00547 real_type t = real_type (); 00548 while (-- size >= 0) { 00549 real_type u (type_traits<value_type>::norm_inf (*it)); 00550 if (u > t) 00551 t = u; 00552 ++ it; 00553 } 00554 return t; 00555 } 00556 // Sparse case 00557 template<class I> 00558 static BOOST_UBLAS_INLINE 00559 result_type apply (I it, const I &it_end) { 00560 real_type t = real_type (); 00561 while (it != it_end) { 00562 real_type u (type_traits<value_type>::norm_inf (*it)); 00563 if (u > t) 00564 t = u; 00565 ++ it; 00566 } 00567 return t; 00568 } 00569 }; 00570 00571 // Unary returning index 00572 template<class V> 00573 struct vector_scalar_index_unary_functor { 00574 typedef typename V::value_type value_type; 00575 typedef typename type_traits<value_type>::real_type real_type; 00576 typedef typename V::size_type result_type; 00577 }; 00578 00579 template<class V> 00580 struct vector_index_norm_inf: 00581 public vector_scalar_index_unary_functor<V> { 00582 typedef typename vector_scalar_index_unary_functor<V>::value_type value_type; 00583 typedef typename vector_scalar_index_unary_functor<V>::real_type real_type; 00584 typedef typename vector_scalar_index_unary_functor<V>::result_type result_type; 00585 00586 template<class E> 00587 static BOOST_UBLAS_INLINE 00588 result_type apply (const vector_expression<E> &e) { 00589 // ISSUE For CBLAS compatibility return 0 index in empty case 00590 result_type i_norm_inf (0); 00591 real_type t = real_type (); 00592 typedef typename E::size_type vector_size_type; 00593 vector_size_type size (e ().size ()); 00594 for (vector_size_type i = 0; i < size; ++ i) { 00595 real_type u (type_traits<value_type>::norm_inf (e () (i))); 00596 if (u > t) { 00597 i_norm_inf = i; 00598 t = u; 00599 } 00600 } 00601 return i_norm_inf; 00602 } 00603 // Dense case 00604 template<class D, class I> 00605 static BOOST_UBLAS_INLINE 00606 result_type apply (D size, I it) { 00607 // ISSUE For CBLAS compatibility return 0 index in empty case 00608 result_type i_norm_inf (0); 00609 real_type t = real_type (); 00610 while (-- size >= 0) { 00611 real_type u (type_traits<value_type>::norm_inf (*it)); 00612 if (u > t) { 00613 i_norm_inf = it.index (); 00614 t = u; 00615 } 00616 ++ it; 00617 } 00618 return i_norm_inf; 00619 } 00620 // Sparse case 00621 template<class I> 00622 static BOOST_UBLAS_INLINE 00623 result_type apply (I it, const I &it_end) { 00624 // ISSUE For CBLAS compatibility return 0 index in empty case 00625 result_type i_norm_inf (0); 00626 real_type t = real_type (); 00627 while (it != it_end) { 00628 real_type u (type_traits<value_type>::norm_inf (*it)); 00629 if (u > t) { 00630 i_norm_inf = it.index (); 00631 t = u; 00632 } 00633 ++ it; 00634 } 00635 return i_norm_inf; 00636 } 00637 }; 00638 00639 // Binary returning scalar 00640 template<class V1, class V2, class TV> 00641 struct vector_scalar_binary_functor { 00642 typedef TV value_type; 00643 typedef TV result_type; 00644 }; 00645 00646 template<class V1, class V2, class TV> 00647 struct vector_inner_prod: 00648 public vector_scalar_binary_functor<V1, V2, TV> { 00649 typedef typename vector_scalar_binary_functor<V1, V2, TV>::value_type value_type; 00650 typedef typename vector_scalar_binary_functor<V1, V2, TV>::result_type result_type; 00651 00652 template<class C1, class C2> 00653 static BOOST_UBLAS_INLINE 00654 result_type apply (const vector_container<C1> &c1, 00655 const vector_container<C2> &c2) { 00656 #ifdef BOOST_UBLAS_USE_SIMD 00657 using namespace raw; 00658 typedef typename C1::size_type vector_size_type; 00659 vector_size_type size (BOOST_UBLAS_SAME (c1 ().size (), c2 ().size ())); 00660 const typename V1::value_type *data1 = data_const (c1 ()); 00661 const typename V1::value_type *data2 = data_const (c2 ()); 00662 vector_size_type s1 = stride (c1 ()); 00663 vector_size_type s2 = stride (c2 ()); 00664 result_type t = result_type (0); 00665 if (s1 == 1 && s2 == 1) { 00666 for (vector_size_type i = 0; i < size; ++ i) 00667 t += data1 [i] * data2 [i]; 00668 } else if (s2 == 1) { 00669 for (vector_size_type i = 0, i1 = 0; i < size; ++ i, i1 += s1) 00670 t += data1 [i1] * data2 [i]; 00671 } else if (s1 == 1) { 00672 for (vector_size_type i = 0, i2 = 0; i < size; ++ i, i2 += s2) 00673 t += data1 [i] * data2 [i2]; 00674 } else { 00675 for (vector_size_type i = 0, i1 = 0, i2 = 0; i < size; ++ i, i1 += s1, i2 += s2) 00676 t += data1 [i1] * data2 [i2]; 00677 } 00678 return t; 00679 #elif defined(BOOST_UBLAS_HAVE_BINDINGS) 00680 return boost::numeric::bindings::atlas::dot (c1 (), c2 ()); 00681 #else 00682 return apply (static_cast<const vector_expression<C1> > (c1), static_cast<const vector_expression<C2> > (c2)); 00683 #endif 00684 } 00685 template<class E1, class E2> 00686 static BOOST_UBLAS_INLINE 00687 result_type apply (const vector_expression<E1> &e1, 00688 const vector_expression<E2> &e2) { 00689 typedef typename E1::size_type vector_size_type; 00690 vector_size_type size (BOOST_UBLAS_SAME (e1 ().size (), e2 ().size ())); 00691 result_type t = result_type (0); 00692 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE 00693 for (vector_size_type i = 0; i < size; ++ i) 00694 t += e1 () (i) * e2 () (i); 00695 #else 00696 vector_size_type i (0); 00697 DD (size, 4, r, (t += e1 () (i) * e2 () (i), ++ i)); 00698 #endif 00699 return t; 00700 } 00701 // Dense case 00702 template<class D, class I1, class I2> 00703 static BOOST_UBLAS_INLINE 00704 result_type apply (D size, I1 it1, I2 it2) { 00705 result_type t = result_type (0); 00706 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE 00707 while (-- size >= 0) 00708 t += *it1 * *it2, ++ it1, ++ it2; 00709 #else 00710 DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2)); 00711 #endif 00712 return t; 00713 } 00714 // Packed case 00715 template<class I1, class I2> 00716 static BOOST_UBLAS_INLINE 00717 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) { 00718 result_type t = result_type (0); 00719 typedef typename I1::difference_type vector_difference_type; 00720 vector_difference_type it1_size (it1_end - it1); 00721 vector_difference_type it2_size (it2_end - it2); 00722 vector_difference_type diff (0); 00723 if (it1_size > 0 && it2_size > 0) 00724 diff = it2.index () - it1.index (); 00725 if (diff != 0) { 00726 vector_difference_type size = (std::min) (diff, it1_size); 00727 if (size > 0) { 00728 it1 += size; 00729 it1_size -= size; 00730 diff -= size; 00731 } 00732 size = (std::min) (- diff, it2_size); 00733 if (size > 0) { 00734 it2 += size; 00735 it2_size -= size; 00736 diff += size; 00737 } 00738 } 00739 vector_difference_type size ((std::min) (it1_size, it2_size)); 00740 while (-- size >= 0) 00741 t += *it1 * *it2, ++ it1, ++ it2; 00742 return t; 00743 } 00744 // Sparse case 00745 template<class I1, class I2> 00746 static BOOST_UBLAS_INLINE 00747 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) { 00748 result_type t = result_type (0); 00749 if (it1 != it1_end && it2 != it2_end) { 00750 while (true) { 00751 if (it1.index () == it2.index ()) { 00752 t += *it1 * *it2, ++ it1, ++ it2; 00753 if (it1 == it1_end || it2 == it2_end) 00754 break; 00755 } else if (it1.index () < it2.index ()) { 00756 increment (it1, it1_end, it2.index () - it1.index ()); 00757 if (it1 == it1_end) 00758 break; 00759 } else if (it1.index () > it2.index ()) { 00760 increment (it2, it2_end, it1.index () - it2.index ()); 00761 if (it2 == it2_end) 00762 break; 00763 } 00764 } 00765 } 00766 return t; 00767 } 00768 }; 00769 00770 // Matrix functors 00771 00772 // Binary returning vector 00773 template<class M1, class M2, class TV> 00774 struct matrix_vector_binary_functor { 00775 typedef typename M1::size_type size_type; 00776 typedef typename M1::difference_type difference_type; 00777 typedef TV value_type; 00778 typedef TV result_type; 00779 }; 00780 00781 template<class M1, class M2, class TV> 00782 struct matrix_vector_prod1: 00783 public matrix_vector_binary_functor<M1, M2, TV> { 00784 typedef typename matrix_vector_binary_functor<M1, M2, TV>::size_type size_type; 00785 typedef typename matrix_vector_binary_functor<M1, M2, TV>::difference_type difference_type; 00786 typedef typename matrix_vector_binary_functor<M1, M2, TV>::value_type value_type; 00787 typedef typename matrix_vector_binary_functor<M1, M2, TV>::result_type result_type; 00788 00789 template<class C1, class C2> 00790 static BOOST_UBLAS_INLINE 00791 result_type apply (const matrix_container<C1> &c1, 00792 const vector_container<C2> &c2, 00793 size_type i) { 00794 #ifdef BOOST_UBLAS_USE_SIMD 00795 using namespace raw; 00796 size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().size ()); 00797 const typename M1::value_type *data1 = data_const (c1 ()) + i * stride1 (c1 ()); 00798 const typename M2::value_type *data2 = data_const (c2 ()); 00799 size_type s1 = stride2 (c1 ()); 00800 size_type s2 = stride (c2 ()); 00801 result_type t = result_type (0); 00802 if (s1 == 1 && s2 == 1) { 00803 for (size_type j = 0; j < size; ++ j) 00804 t += data1 [j] * data2 [j]; 00805 } else if (s2 == 1) { 00806 for (size_type j = 0, j1 = 0; j < size; ++ j, j1 += s1) 00807 t += data1 [j1] * data2 [j]; 00808 } else if (s1 == 1) { 00809 for (size_type j = 0, j2 = 0; j < size; ++ j, j2 += s2) 00810 t += data1 [j] * data2 [j2]; 00811 } else { 00812 for (size_type j = 0, j1 = 0, j2 = 0; j < size; ++ j, j1 += s1, j2 += s2) 00813 t += data1 [j1] * data2 [j2]; 00814 } 00815 return t; 00816 #elif defined(BOOST_UBLAS_HAVE_BINDINGS) 00817 return boost::numeric::bindings::atlas::dot (c1 ().row (i), c2 ()); 00818 #else 00819 return apply (static_cast<const matrix_expression<C1> > (c1), static_cast<const vector_expression<C2> > (c2, i)); 00820 #endif 00821 } 00822 template<class E1, class E2> 00823 static BOOST_UBLAS_INLINE 00824 result_type apply (const matrix_expression<E1> &e1, 00825 const vector_expression<E2> &e2, 00826 size_type i) { 00827 size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size ()); 00828 result_type t = result_type (0); 00829 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE 00830 for (size_type j = 0; j < size; ++ j) 00831 t += e1 () (i, j) * e2 () (j); 00832 #else 00833 size_type j (0); 00834 DD (size, 4, r, (t += e1 () (i, j) * e2 () (j), ++ j)); 00835 #endif 00836 return t; 00837 } 00838 // Dense case 00839 template<class I1, class I2> 00840 static BOOST_UBLAS_INLINE 00841 result_type apply (difference_type size, I1 it1, I2 it2) { 00842 result_type t = result_type (0); 00843 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE 00844 while (-- size >= 0) 00845 t += *it1 * *it2, ++ it1, ++ it2; 00846 #else 00847 DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2)); 00848 #endif 00849 return t; 00850 } 00851 // Packed case 00852 template<class I1, class I2> 00853 static BOOST_UBLAS_INLINE 00854 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) { 00855 result_type t = result_type (0); 00856 difference_type it1_size (it1_end - it1); 00857 difference_type it2_size (it2_end - it2); 00858 difference_type diff (0); 00859 if (it1_size > 0 && it2_size > 0) 00860 diff = it2.index () - it1.index2 (); 00861 if (diff != 0) { 00862 difference_type size = (std::min) (diff, it1_size); 00863 if (size > 0) { 00864 it1 += size; 00865 it1_size -= size; 00866 diff -= size; 00867 } 00868 size = (std::min) (- diff, it2_size); 00869 if (size > 0) { 00870 it2 += size; 00871 it2_size -= size; 00872 diff += size; 00873 } 00874 } 00875 difference_type size ((std::min) (it1_size, it2_size)); 00876 while (-- size >= 0) 00877 t += *it1 * *it2, ++ it1, ++ it2; 00878 return t; 00879 } 00880 // Sparse case 00881 template<class I1, class I2> 00882 static BOOST_UBLAS_INLINE 00883 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, 00884 sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) { 00885 result_type t = result_type (0); 00886 if (it1 != it1_end && it2 != it2_end) { 00887 size_type it1_index = it1.index2 (), it2_index = it2.index (); 00888 while (true) { 00889 difference_type compare = it1_index - it2_index; 00890 if (compare == 0) { 00891 t += *it1 * *it2, ++ it1, ++ it2; 00892 if (it1 != it1_end && it2 != it2_end) { 00893 it1_index = it1.index2 (); 00894 it2_index = it2.index (); 00895 } else 00896 break; 00897 } else if (compare < 0) { 00898 increment (it1, it1_end, - compare); 00899 if (it1 != it1_end) 00900 it1_index = it1.index2 (); 00901 else 00902 break; 00903 } else if (compare > 0) { 00904 increment (it2, it2_end, compare); 00905 if (it2 != it2_end) 00906 it2_index = it2.index (); 00907 else 00908 break; 00909 } 00910 } 00911 } 00912 return t; 00913 } 00914 // Sparse packed case 00915 template<class I1, class I2> 00916 static BOOST_UBLAS_INLINE 00917 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */, 00918 sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) { 00919 result_type t = result_type (0); 00920 while (it1 != it1_end) { 00921 t += *it1 * it2 () (it1.index2 ()); 00922 ++ it1; 00923 } 00924 return t; 00925 } 00926 // Packed sparse case 00927 template<class I1, class I2> 00928 static BOOST_UBLAS_INLINE 00929 result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end, 00930 packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) { 00931 result_type t = result_type (0); 00932 while (it2 != it2_end) { 00933 t += it1 () (it1.index1 (), it2.index ()) * *it2; 00934 ++ it2; 00935 } 00936 return t; 00937 } 00938 // Another dispatcher 00939 template<class I1, class I2> 00940 static BOOST_UBLAS_INLINE 00941 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, 00942 sparse_bidirectional_iterator_tag) { 00943 typedef typename I1::iterator_category iterator1_category; 00944 typedef typename I2::iterator_category iterator2_category; 00945 return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ()); 00946 } 00947 }; 00948 00949 template<class M1, class M2, class TV> 00950 struct matrix_vector_prod2: 00951 public matrix_vector_binary_functor<M1, M2, TV> { 00952 typedef typename matrix_vector_binary_functor<M1, M2, TV>::size_type size_type; 00953 typedef typename matrix_vector_binary_functor<M1, M2, TV>::difference_type difference_type; 00954 typedef typename matrix_vector_binary_functor<M1, M2, TV>::value_type value_type; 00955 typedef typename matrix_vector_binary_functor<M1, M2, TV>::result_type result_type; 00956 00957 template<class C1, class C2> 00958 static BOOST_UBLAS_INLINE 00959 result_type apply (const vector_container<C1> &c1, 00960 const matrix_container<C2> &c2, 00961 size_type i) { 00962 #ifdef BOOST_UBLAS_USE_SIMD 00963 using namespace raw; 00964 size_type size = BOOST_UBLAS_SAME (c1 ().size (), c2 ().size1 ()); 00965 const typename M1::value_type *data1 = data_const (c1 ()); 00966 const typename M2::value_type *data2 = data_const (c2 ()) + i * stride2 (c2 ()); 00967 size_type s1 = stride (c1 ()); 00968 size_type s2 = stride1 (c2 ()); 00969 result_type t = result_type (0); 00970 if (s1 == 1 && s2 == 1) { 00971 for (size_type j = 0; j < size; ++ j) 00972 t += data1 [j] * data2 [j]; 00973 } else if (s2 == 1) { 00974 for (size_type j = 0, j1 = 0; j < size; ++ j, j1 += s1) 00975 t += data1 [j1] * data2 [j]; 00976 } else if (s1 == 1) { 00977 for (size_type j = 0, j2 = 0; j < size; ++ j, j2 += s2) 00978 t += data1 [j] * data2 [j2]; 00979 } else { 00980 for (size_type j = 0, j1 = 0, j2 = 0; j < size; ++ j, j1 += s1, j2 += s2) 00981 t += data1 [j1] * data2 [j2]; 00982 } 00983 return t; 00984 #elif defined(BOOST_UBLAS_HAVE_BINDINGS) 00985 return boost::numeric::bindings::atlas::dot (c1 (), c2 ().column (i)); 00986 #else 00987 return apply (static_cast<const vector_expression<C1> > (c1), static_cast<const matrix_expression<C2> > (c2, i)); 00988 #endif 00989 } 00990 template<class E1, class E2> 00991 static BOOST_UBLAS_INLINE 00992 result_type apply (const vector_expression<E1> &e1, 00993 const matrix_expression<E2> &e2, 00994 size_type i) { 00995 size_type size = BOOST_UBLAS_SAME (e1 ().size (), e2 ().size1 ()); 00996 result_type t = result_type (0); 00997 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE 00998 for (size_type j = 0; j < size; ++ j) 00999 t += e1 () (j) * e2 () (j, i); 01000 #else 01001 size_type j (0); 01002 DD (size, 4, r, (t += e1 () (j) * e2 () (j, i), ++ j)); 01003 #endif 01004 return t; 01005 } 01006 // Dense case 01007 template<class I1, class I2> 01008 static BOOST_UBLAS_INLINE 01009 result_type apply (difference_type size, I1 it1, I2 it2) { 01010 result_type t = result_type (0); 01011 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE 01012 while (-- size >= 0) 01013 t += *it1 * *it2, ++ it1, ++ it2; 01014 #else 01015 DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2)); 01016 #endif 01017 return t; 01018 } 01019 // Packed case 01020 template<class I1, class I2> 01021 static BOOST_UBLAS_INLINE 01022 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) { 01023 result_type t = result_type (0); 01024 difference_type it1_size (it1_end - it1); 01025 difference_type it2_size (it2_end - it2); 01026 difference_type diff (0); 01027 if (it1_size > 0 && it2_size > 0) 01028 diff = it2.index1 () - it1.index (); 01029 if (diff != 0) { 01030 difference_type size = (std::min) (diff, it1_size); 01031 if (size > 0) { 01032 it1 += size; 01033 it1_size -= size; 01034 diff -= size; 01035 } 01036 size = (std::min) (- diff, it2_size); 01037 if (size > 0) { 01038 it2 += size; 01039 it2_size -= size; 01040 diff += size; 01041 } 01042 } 01043 difference_type size ((std::min) (it1_size, it2_size)); 01044 while (-- size >= 0) 01045 t += *it1 * *it2, ++ it1, ++ it2; 01046 return t; 01047 } 01048 // Sparse case 01049 template<class I1, class I2> 01050 static BOOST_UBLAS_INLINE 01051 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, 01052 sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) { 01053 result_type t = result_type (0); 01054 if (it1 != it1_end && it2 != it2_end) { 01055 size_type it1_index = it1.index (), it2_index = it2.index1 (); 01056 while (true) { 01057 difference_type compare = it1_index - it2_index; 01058 if (compare == 0) { 01059 t += *it1 * *it2, ++ it1, ++ it2; 01060 if (it1 != it1_end && it2 != it2_end) { 01061 it1_index = it1.index (); 01062 it2_index = it2.index1 (); 01063 } else 01064 break; 01065 } else if (compare < 0) { 01066 increment (it1, it1_end, - compare); 01067 if (it1 != it1_end) 01068 it1_index = it1.index (); 01069 else 01070 break; 01071 } else if (compare > 0) { 01072 increment (it2, it2_end, compare); 01073 if (it2 != it2_end) 01074 it2_index = it2.index1 (); 01075 else 01076 break; 01077 } 01078 } 01079 } 01080 return t; 01081 } 01082 // Packed sparse case 01083 template<class I1, class I2> 01084 static BOOST_UBLAS_INLINE 01085 result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end, 01086 packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) { 01087 result_type t = result_type (0); 01088 while (it2 != it2_end) { 01089 t += it1 () (it2.index1 ()) * *it2; 01090 ++ it2; 01091 } 01092 return t; 01093 } 01094 // Sparse packed case 01095 template<class I1, class I2> 01096 static BOOST_UBLAS_INLINE 01097 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */, 01098 sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) { 01099 result_type t = result_type (0); 01100 while (it1 != it1_end) { 01101 t += *it1 * it2 () (it1.index (), it2.index2 ()); 01102 ++ it1; 01103 } 01104 return t; 01105 } 01106 // Another dispatcher 01107 template<class I1, class I2> 01108 static BOOST_UBLAS_INLINE 01109 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, 01110 sparse_bidirectional_iterator_tag) { 01111 typedef typename I1::iterator_category iterator1_category; 01112 typedef typename I2::iterator_category iterator2_category; 01113 return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ()); 01114 } 01115 }; 01116 01117 // Binary returning matrix 01118 template<class M1, class M2, class TV> 01119 struct matrix_matrix_binary_functor { 01120 typedef typename M1::size_type size_type; 01121 typedef typename M1::difference_type difference_type; 01122 typedef TV value_type; 01123 typedef TV result_type; 01124 }; 01125 01126 template<class M1, class M2, class TV> 01127 struct matrix_matrix_prod: 01128 public matrix_matrix_binary_functor<M1, M2, TV> { 01129 typedef typename matrix_matrix_binary_functor<M1, M2, TV>::size_type size_type; 01130 typedef typename matrix_matrix_binary_functor<M1, M2, TV>::difference_type difference_type; 01131 typedef typename matrix_matrix_binary_functor<M1, M2, TV>::value_type value_type; 01132 typedef typename matrix_matrix_binary_functor<M1, M2, TV>::result_type result_type; 01133 01134 template<class C1, class C2> 01135 static BOOST_UBLAS_INLINE 01136 result_type apply (const matrix_container<C1> &c1, 01137 const matrix_container<C2> &c2, 01138 size_type i, size_type j) { 01139 #ifdef BOOST_UBLAS_USE_SIMD 01140 using namespace raw; 01141 size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().sizc1 ()); 01142 const typename M1::value_type *data1 = data_const (c1 ()) + i * stride1 (c1 ()); 01143 const typename M2::value_type *data2 = data_const (c2 ()) + j * stride2 (c2 ()); 01144 size_type s1 = stride2 (c1 ()); 01145 size_type s2 = stride1 (c2 ()); 01146 result_type t = result_type (0); 01147 if (s1 == 1 && s2 == 1) { 01148 for (size_type k = 0; k < size; ++ k) 01149 t += data1 [k] * data2 [k]; 01150 } else if (s2 == 1) { 01151 for (size_type k = 0, k1 = 0; k < size; ++ k, k1 += s1) 01152 t += data1 [k1] * data2 [k]; 01153 } else if (s1 == 1) { 01154 for (size_type k = 0, k2 = 0; k < size; ++ k, k2 += s2) 01155 t += data1 [k] * data2 [k2]; 01156 } else { 01157 for (size_type k = 0, k1 = 0, k2 = 0; k < size; ++ k, k1 += s1, k2 += s2) 01158 t += data1 [k1] * data2 [k2]; 01159 } 01160 return t; 01161 #elif defined(BOOST_UBLAS_HAVE_BINDINGS) 01162 return boost::numeric::bindings::atlas::dot (c1 ().row (i), c2 ().column (j)); 01163 #else 01164 return apply (static_cast<const matrix_expression<C1> > (c1), static_cast<const matrix_expression<C2> > (c2, i)); 01165 #endif 01166 } 01167 template<class E1, class E2> 01168 static BOOST_UBLAS_INLINE 01169 result_type apply (const matrix_expression<E1> &e1, 01170 const matrix_expression<E2> &e2, 01171 size_type i, size_type j) { 01172 size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ()); 01173 result_type t = result_type (0); 01174 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE 01175 for (size_type k = 0; k < size; ++ k) 01176 t += e1 () (i, k) * e2 () (k, j); 01177 #else 01178 size_type k (0); 01179 DD (size, 4, r, (t += e1 () (i, k) * e2 () (k, j), ++ k)); 01180 #endif 01181 return t; 01182 } 01183 // Dense case 01184 template<class I1, class I2> 01185 static BOOST_UBLAS_INLINE 01186 result_type apply (difference_type size, I1 it1, I2 it2) { 01187 result_type t = result_type (0); 01188 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE 01189 while (-- size >= 0) 01190 t += *it1 * *it2, ++ it1, ++ it2; 01191 #else 01192 DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2)); 01193 #endif 01194 return t; 01195 } 01196 // Packed case 01197 template<class I1, class I2> 01198 static BOOST_UBLAS_INLINE 01199 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, packed_random_access_iterator_tag) { 01200 result_type t = result_type (0); 01201 difference_type it1_size (it1_end - it1); 01202 difference_type it2_size (it2_end - it2); 01203 difference_type diff (0); 01204 if (it1_size > 0 && it2_size > 0) 01205 diff = it2.index1 () - it1.index2 (); 01206 if (diff != 0) { 01207 difference_type size = (std::min) (diff, it1_size); 01208 if (size > 0) { 01209 it1 += size; 01210 it1_size -= size; 01211 diff -= size; 01212 } 01213 size = (std::min) (- diff, it2_size); 01214 if (size > 0) { 01215 it2 += size; 01216 it2_size -= size; 01217 diff += size; 01218 } 01219 } 01220 difference_type size ((std::min) (it1_size, it2_size)); 01221 while (-- size >= 0) 01222 t += *it1 * *it2, ++ it1, ++ it2; 01223 return t; 01224 } 01225 // Sparse case 01226 template<class I1, class I2> 01227 static BOOST_UBLAS_INLINE 01228 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) { 01229 result_type t = result_type (0); 01230 if (it1 != it1_end && it2 != it2_end) { 01231 size_type it1_index = it1.index2 (), it2_index = it2.index1 (); 01232 while (true) { 01233 difference_type compare = difference_type (it1_index - it2_index); 01234 if (compare == 0) { 01235 t += *it1 * *it2, ++ it1, ++ it2; 01236 if (it1 != it1_end && it2 != it2_end) { 01237 it1_index = it1.index2 (); 01238 it2_index = it2.index1 (); 01239 } else 01240 break; 01241 } else if (compare < 0) { 01242 increment (it1, it1_end, - compare); 01243 if (it1 != it1_end) 01244 it1_index = it1.index2 (); 01245 else 01246 break; 01247 } else if (compare > 0) { 01248 increment (it2, it2_end, compare); 01249 if (it2 != it2_end) 01250 it2_index = it2.index1 (); 01251 else 01252 break; 01253 } 01254 } 01255 } 01256 return t; 01257 } 01258 }; 01259 01260 // Unary returning scalar norm 01261 template<class M> 01262 struct matrix_scalar_real_unary_functor { 01263 typedef typename M::value_type value_type; 01264 typedef typename type_traits<value_type>::real_type real_type; 01265 typedef real_type result_type; 01266 }; 01267 01268 template<class M> 01269 struct matrix_norm_1: 01270 public matrix_scalar_real_unary_functor<M> { 01271 typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type; 01272 typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type; 01273 typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type; 01274 01275 template<class E> 01276 static BOOST_UBLAS_INLINE 01277 result_type apply (const matrix_expression<E> &e) { 01278 real_type t = real_type (); 01279 typedef typename E::size_type matrix_size_type; 01280 matrix_size_type size2 (e ().size2 ()); 01281 for (matrix_size_type j = 0; j < size2; ++ j) { 01282 real_type u = real_type (); 01283 matrix_size_type size1 (e ().size1 ()); 01284 for (matrix_size_type i = 0; i < size1; ++ i) { 01285 real_type v (type_traits<value_type>::norm_1 (e () (i, j))); 01286 u += v; 01287 } 01288 if (u > t) 01289 t = u; 01290 } 01291 return t; 01292 } 01293 }; 01294 01295 template<class M> 01296 struct matrix_norm_frobenius: 01297 public matrix_scalar_real_unary_functor<M> { 01298 typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type; 01299 typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type; 01300 typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type; 01301 01302 template<class E> 01303 static BOOST_UBLAS_INLINE 01304 result_type apply (const matrix_expression<E> &e) { 01305 real_type t = real_type (); 01306 typedef typename E::size_type matrix_size_type; 01307 matrix_size_type size1 (e ().size1 ()); 01308 for (matrix_size_type i = 0; i < size1; ++ i) { 01309 matrix_size_type size2 (e ().size2 ()); 01310 for (matrix_size_type j = 0; j < size2; ++ j) { 01311 real_type u (type_traits<value_type>::norm_2 (e () (i, j))); 01312 t += u * u; 01313 } 01314 } 01315 return type_traits<real_type>::type_sqrt (t); 01316 } 01317 }; 01318 01319 template<class M> 01320 struct matrix_norm_inf: 01321 public matrix_scalar_real_unary_functor<M> { 01322 typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type; 01323 typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type; 01324 typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type; 01325 01326 template<class E> 01327 static BOOST_UBLAS_INLINE 01328 result_type apply (const matrix_expression<E> &e) { 01329 real_type t = real_type (); 01330 typedef typename E::size_type matrix_size_type; 01331 matrix_size_type size1 (e ().size1 ()); 01332 for (matrix_size_type i = 0; i < size1; ++ i) { 01333 real_type u = real_type (); 01334 matrix_size_type size2 (e ().size2 ()); 01335 for (matrix_size_type j = 0; j < size2; ++ j) { 01336 real_type v (type_traits<value_type>::norm_inf (e () (i, j))); 01337 u += v; 01338 } 01339 if (u > t) 01340 t = u; 01341 } 01342 return t; 01343 } 01344 }; 01345 01346 // forward declaration 01347 template <class Z, class D> struct basic_column_major; 01348 01349 // This functor defines storage layout and it's properties 01350 // matrix (i,j) -> storage [i * size_i + j] 01351 template <class Z, class D> 01352 struct basic_row_major { 01353 typedef Z size_type; 01354 typedef D difference_type; 01355 typedef row_major_tag orientation_category; 01356 typedef basic_column_major<Z,D> transposed_layout; 01357 01358 static 01359 BOOST_UBLAS_INLINE 01360 size_type storage_size (size_type size_i, size_type size_j) { 01361 // Guard against size_type overflow 01362 BOOST_UBLAS_CHECK (size_j == 0 || size_i <= (std::numeric_limits<size_type>::max) () / size_j, bad_size ()); 01363 return size_i * size_j; 01364 } 01365 01366 // Indexing conversion to storage element 01367 static 01368 BOOST_UBLAS_INLINE 01369 size_type element (size_type i, size_type size_i, size_type j, size_type size_j) { 01370 BOOST_UBLAS_CHECK (i < size_i, bad_index ()); 01371 BOOST_UBLAS_CHECK (j < size_j, bad_index ()); 01372 detail::ignore_unused_variable_warning(size_i); 01373 // Guard against size_type overflow 01374 BOOST_UBLAS_CHECK (i <= ((std::numeric_limits<size_type>::max) () - j) / size_j, bad_index ()); 01375 return i * size_j + j; 01376 } 01377 static 01378 BOOST_UBLAS_INLINE 01379 size_type address (size_type i, size_type size_i, size_type j, size_type size_j) { 01380 BOOST_UBLAS_CHECK (i <= size_i, bad_index ()); 01381 BOOST_UBLAS_CHECK (j <= size_j, bad_index ()); 01382 // Guard against size_type overflow - address may be size_j past end of storage 01383 BOOST_UBLAS_CHECK (size_j == 0 || i <= ((std::numeric_limits<size_type>::max) () - j) / size_j, bad_index ()); 01384 detail::ignore_unused_variable_warning(size_i); 01385 return i * size_j + j; 01386 } 01387 01388 // Storage element to index conversion 01389 static 01390 BOOST_UBLAS_INLINE 01391 difference_type distance_i (difference_type k, size_type /* size_i */, size_type size_j) { 01392 return size_j != 0 ? k / size_j : 0; 01393 } 01394 static 01395 BOOST_UBLAS_INLINE 01396 difference_type distance_j (difference_type k, size_type /* size_i */, size_type /* size_j */) { 01397 return k; 01398 } 01399 static 01400 BOOST_UBLAS_INLINE 01401 size_type index_i (difference_type k, size_type /* size_i */, size_type size_j) { 01402 return size_j != 0 ? k / size_j : 0; 01403 } 01404 static 01405 BOOST_UBLAS_INLINE 01406 size_type index_j (difference_type k, size_type /* size_i */, size_type size_j) { 01407 return size_j != 0 ? k % size_j : 0; 01408 } 01409 static 01410 BOOST_UBLAS_INLINE 01411 bool fast_i () { 01412 return false; 01413 } 01414 static 01415 BOOST_UBLAS_INLINE 01416 bool fast_j () { 01417 return true; 01418 } 01419 01420 // Iterating storage elements 01421 template<class I> 01422 static 01423 BOOST_UBLAS_INLINE 01424 void increment_i (I &it, size_type /* size_i */, size_type size_j) { 01425 it += size_j; 01426 } 01427 template<class I> 01428 static 01429 BOOST_UBLAS_INLINE 01430 void increment_i (I &it, difference_type n, size_type /* size_i */, size_type size_j) { 01431 it += n * size_j; 01432 } 01433 template<class I> 01434 static 01435 BOOST_UBLAS_INLINE 01436 void decrement_i (I &it, size_type /* size_i */, size_type size_j) { 01437 it -= size_j; 01438 } 01439 template<class I> 01440 static 01441 BOOST_UBLAS_INLINE 01442 void decrement_i (I &it, difference_type n, size_type /* size_i */, size_type size_j) { 01443 it -= n * size_j; 01444 } 01445 template<class I> 01446 static 01447 BOOST_UBLAS_INLINE 01448 void increment_j (I &it, size_type /* size_i */, size_type /* size_j */) { 01449 ++ it; 01450 } 01451 template<class I> 01452 static 01453 BOOST_UBLAS_INLINE 01454 void increment_j (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) { 01455 it += n; 01456 } 01457 template<class I> 01458 static 01459 BOOST_UBLAS_INLINE 01460 void decrement_j (I &it, size_type /* size_i */, size_type /* size_j */) { 01461 -- it; 01462 } 01463 template<class I> 01464 static 01465 BOOST_UBLAS_INLINE 01466 void decrement_j (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) { 01467 it -= n; 01468 } 01469 01470 // Triangular access 01471 static 01472 BOOST_UBLAS_INLINE 01473 size_type triangular_size (size_type size_i, size_type size_j) { 01474 size_type size = (std::max) (size_i, size_j); 01475 // Guard against size_type overflow - simplified 01476 BOOST_UBLAS_CHECK (size == 0 || size / 2 < (std::numeric_limits<size_type>::max) () / size /* +1/2 */, bad_size ()); 01477 return ((size + 1) * size) / 2; 01478 } 01479 static 01480 BOOST_UBLAS_INLINE 01481 size_type lower_element (size_type i, size_type size_i, size_type j, size_type size_j) { 01482 BOOST_UBLAS_CHECK (i < size_i, bad_index ()); 01483 BOOST_UBLAS_CHECK (j < size_j, bad_index ()); 01484 BOOST_UBLAS_CHECK (i >= j, bad_index ()); 01485 detail::ignore_unused_variable_warning(size_i); 01486 detail::ignore_unused_variable_warning(size_j); 01487 // FIXME size_type overflow 01488 // sigma_i (i + 1) = (i + 1) * i / 2 01489 // i = 0 1 2 3, sigma = 0 1 3 6 01490 return ((i + 1) * i) / 2 + j; 01491 } 01492 static 01493 BOOST_UBLAS_INLINE 01494 size_type upper_element (size_type i, size_type size_i, size_type j, size_type size_j) { 01495 BOOST_UBLAS_CHECK (i < size_i, bad_index ()); 01496 BOOST_UBLAS_CHECK (j < size_j, bad_index ()); 01497 BOOST_UBLAS_CHECK (i <= j, bad_index ()); 01498 // FIXME size_type overflow 01499 // sigma_i (size - i) = size * i - i * (i - 1) / 2 01500 // i = 0 1 2 3, sigma = 0 4 7 9 01501 return (i * (2 * (std::max) (size_i, size_j) - i + 1)) / 2 + j - i; 01502 } 01503 01504 // Major and minor indices 01505 static 01506 BOOST_UBLAS_INLINE 01507 size_type index_M (size_type index1, size_type /* index2 */) { 01508 return index1; 01509 } 01510 static 01511 BOOST_UBLAS_INLINE 01512 size_type index_m (size_type /* index1 */, size_type index2) { 01513 return index2; 01514 } 01515 static 01516 BOOST_UBLAS_INLINE 01517 size_type size_M (size_type size_i, size_type /* size_j */) { 01518 return size_i; 01519 } 01520 static 01521 BOOST_UBLAS_INLINE 01522 size_type size_m (size_type /* size_i */, size_type size_j) { 01523 return size_j; 01524 } 01525 }; 01526 01527 // This functor defines storage layout and it's properties 01528 // matrix (i,j) -> storage [i + j * size_i] 01529 template <class Z, class D> 01530 struct basic_column_major { 01531 typedef Z size_type; 01532 typedef D difference_type; 01533 typedef column_major_tag orientation_category; 01534 typedef basic_row_major<Z,D> transposed_layout; 01535 01536 static 01537 BOOST_UBLAS_INLINE 01538 size_type storage_size (size_type size_i, size_type size_j) { 01539 // Guard against size_type overflow 01540 BOOST_UBLAS_CHECK (size_i == 0 || size_j <= (std::numeric_limits<size_type>::max) () / size_i, bad_size ()); 01541 return size_i * size_j; 01542 } 01543 01544 // Indexing conversion to storage element 01545 static 01546 BOOST_UBLAS_INLINE 01547 size_type element (size_type i, size_type size_i, size_type j, size_type size_j) { 01548 BOOST_UBLAS_CHECK (i < size_i, bad_index ()); 01549 BOOST_UBLAS_CHECK (j < size_j, bad_index ()); 01550 detail::ignore_unused_variable_warning(size_j); 01551 // Guard against size_type overflow 01552 BOOST_UBLAS_CHECK (j <= ((std::numeric_limits<size_type>::max) () - i) / size_i, bad_index ()); 01553 return i + j * size_i; 01554 } 01555 static 01556 BOOST_UBLAS_INLINE 01557 size_type address (size_type i, size_type size_i, size_type j, size_type size_j) { 01558 BOOST_UBLAS_CHECK (i <= size_i, bad_index ()); 01559 BOOST_UBLAS_CHECK (j <= size_j, bad_index ()); 01560 detail::ignore_unused_variable_warning(size_j); 01561 // Guard against size_type overflow - address may be size_i past end of storage 01562 BOOST_UBLAS_CHECK (size_i == 0 || j <= ((std::numeric_limits<size_type>::max) () - i) / size_i, bad_index ()); 01563 return i + j * size_i; 01564 } 01565 01566 // Storage element to index conversion 01567 static 01568 BOOST_UBLAS_INLINE 01569 difference_type distance_i (difference_type k, size_type /* size_i */, size_type /* size_j */) { 01570 return k; 01571 } 01572 static 01573 BOOST_UBLAS_INLINE 01574 difference_type distance_j (difference_type k, size_type size_i, size_type /* size_j */) { 01575 return size_i != 0 ? k / size_i : 0; 01576 } 01577 static 01578 BOOST_UBLAS_INLINE 01579 size_type index_i (difference_type k, size_type size_i, size_type /* size_j */) { 01580 return size_i != 0 ? k % size_i : 0; 01581 } 01582 static 01583 BOOST_UBLAS_INLINE 01584 size_type index_j (difference_type k, size_type size_i, size_type /* size_j */) { 01585 return size_i != 0 ? k / size_i : 0; 01586 } 01587 static 01588 BOOST_UBLAS_INLINE 01589 bool fast_i () { 01590 return true; 01591 } 01592 static 01593 BOOST_UBLAS_INLINE 01594 bool fast_j () { 01595 return false; 01596 } 01597 01598 // Iterating 01599 template<class I> 01600 static 01601 BOOST_UBLAS_INLINE 01602 void increment_i (I &it, size_type /* size_i */, size_type /* size_j */) { 01603 ++ it; 01604 } 01605 template<class I> 01606 static 01607 BOOST_UBLAS_INLINE 01608 void increment_i (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) { 01609 it += n; 01610 } 01611 template<class I> 01612 static 01613 BOOST_UBLAS_INLINE 01614 void decrement_i (I &it, size_type /* size_i */, size_type /* size_j */) { 01615 -- it; 01616 } 01617 template<class I> 01618 static 01619 BOOST_UBLAS_INLINE 01620 void decrement_i (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) { 01621 it -= n; 01622 } 01623 template<class I> 01624 static 01625 BOOST_UBLAS_INLINE 01626 void increment_j (I &it, size_type size_i, size_type /* size_j */) { 01627 it += size_i; 01628 } 01629 template<class I> 01630 static 01631 BOOST_UBLAS_INLINE 01632 void increment_j (I &it, difference_type n, size_type size_i, size_type /* size_j */) { 01633 it += n * size_i; 01634 } 01635 template<class I> 01636 static 01637 BOOST_UBLAS_INLINE 01638 void decrement_j (I &it, size_type size_i, size_type /* size_j */) { 01639 it -= size_i; 01640 } 01641 template<class I> 01642 static 01643 BOOST_UBLAS_INLINE 01644 void decrement_j (I &it, difference_type n, size_type size_i, size_type /* size_j */) { 01645 it -= n* size_i; 01646 } 01647 01648 // Triangular access 01649 static 01650 BOOST_UBLAS_INLINE 01651 size_type triangular_size (size_type size_i, size_type size_j) { 01652 size_type size = (std::max) (size_i, size_j); 01653 // Guard against size_type overflow - simplified 01654 BOOST_UBLAS_CHECK (size == 0 || size / 2 < (std::numeric_limits<size_type>::max) () / size /* +1/2 */, bad_size ()); 01655 return ((size + 1) * size) / 2; 01656 } 01657 static 01658 BOOST_UBLAS_INLINE 01659 size_type lower_element (size_type i, size_type size_i, size_type j, size_type size_j) { 01660 BOOST_UBLAS_CHECK (i < size_i, bad_index ()); 01661 BOOST_UBLAS_CHECK (j < size_j, bad_index ()); 01662 BOOST_UBLAS_CHECK (i >= j, bad_index ()); 01663 // FIXME size_type overflow 01664 // sigma_j (size - j) = size * j - j * (j - 1) / 2 01665 // j = 0 1 2 3, sigma = 0 4 7 9 01666 return i - j + (j * (2 * (std::max) (size_i, size_j) - j + 1)) / 2; 01667 } 01668 static 01669 BOOST_UBLAS_INLINE 01670 size_type upper_element (size_type i, size_type size_i, size_type j, size_type size_j) { 01671 BOOST_UBLAS_CHECK (i < size_i, bad_index ()); 01672 BOOST_UBLAS_CHECK (j < size_j, bad_index ()); 01673 BOOST_UBLAS_CHECK (i <= j, bad_index ()); 01674 // FIXME size_type overflow 01675 // sigma_j (j + 1) = (j + 1) * j / 2 01676 // j = 0 1 2 3, sigma = 0 1 3 6 01677 return i + ((j + 1) * j) / 2; 01678 } 01679 01680 // Major and minor indices 01681 static 01682 BOOST_UBLAS_INLINE 01683 size_type index_M (size_type /* index1 */, size_type index2) { 01684 return index2; 01685 } 01686 static 01687 BOOST_UBLAS_INLINE 01688 size_type index_m (size_type index1, size_type /* index2 */) { 01689 return index1; 01690 } 01691 static 01692 BOOST_UBLAS_INLINE 01693 size_type size_M (size_type /* size_i */, size_type size_j) { 01694 return size_j; 01695 } 01696 static 01697 BOOST_UBLAS_INLINE 01698 size_type size_m (size_type size_i, size_type /* size_j */) { 01699 return size_i; 01700 } 01701 }; 01702 01703 01704 template <class Z> 01705 struct basic_full { 01706 typedef Z size_type; 01707 01708 template<class L> 01709 static 01710 BOOST_UBLAS_INLINE 01711 size_type packed_size (L, size_type size_i, size_type size_j) { 01712 return L::storage_size (size_i, size_j); 01713 } 01714 01715 static 01716 BOOST_UBLAS_INLINE 01717 bool zero (size_type /* i */, size_type /* j */) { 01718 return false; 01719 } 01720 static 01721 BOOST_UBLAS_INLINE 01722 bool one (size_type /* i */, size_type /* j */) { 01723 return false; 01724 } 01725 static 01726 BOOST_UBLAS_INLINE 01727 bool other (size_type /* i */, size_type /* j */) { 01728 return true; 01729 } 01730 // FIXME: this should not be used at all 01731 static 01732 BOOST_UBLAS_INLINE 01733 size_type restrict1 (size_type i, size_type /* j */) { 01734 return i; 01735 } 01736 static 01737 BOOST_UBLAS_INLINE 01738 size_type restrict2 (size_type /* i */, size_type j) { 01739 return j; 01740 } 01741 static 01742 BOOST_UBLAS_INLINE 01743 size_type mutable_restrict1 (size_type i, size_type /* j */) { 01744 return i; 01745 } 01746 static 01747 BOOST_UBLAS_INLINE 01748 size_type mutable_restrict2 (size_type /* i */, size_type j) { 01749 return j; 01750 } 01751 }; 01752 01753 namespace detail { 01754 template < class L > 01755 struct transposed_structure { 01756 typedef typename L::size_type size_type; 01757 01758 template<class LAYOUT> 01759 static 01760 BOOST_UBLAS_INLINE 01761 size_type packed_size (LAYOUT l, size_type size_i, size_type size_j) { 01762 return L::packed_size(l, size_j, size_i); 01763 } 01764 01765 static 01766 BOOST_UBLAS_INLINE 01767 bool zero (size_type i, size_type j) { 01768 return L::zero(j, i); 01769 } 01770 static 01771 BOOST_UBLAS_INLINE 01772 bool one (size_type i, size_type j) { 01773 return L::one(j, i); 01774 } 01775 static 01776 BOOST_UBLAS_INLINE 01777 bool other (size_type i, size_type j) { 01778 return L::other(j, i); 01779 } 01780 template<class LAYOUT> 01781 static 01782 BOOST_UBLAS_INLINE 01783 size_type element (LAYOUT /* l */, size_type i, size_type size_i, size_type j, size_type size_j) { 01784 return L::element(typename LAYOUT::transposed_layout(), j, size_j, i, size_i); 01785 } 01786 01787 static 01788 BOOST_UBLAS_INLINE 01789 size_type restrict1 (size_type i, size_type j, size_type size1, size_type size2) { 01790 return L::restrict2(j, i, size2, size1); 01791 } 01792 static 01793 BOOST_UBLAS_INLINE 01794 size_type restrict2 (size_type i, size_type j, size_type size1, size_type size2) { 01795 return L::restrict1(j, i, size2, size1); 01796 } 01797 static 01798 BOOST_UBLAS_INLINE 01799 size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type size2) { 01800 return L::mutable_restrict2(j, i, size2, size1); 01801 } 01802 static 01803 BOOST_UBLAS_INLINE 01804 size_type mutable_restrict2 (size_type i, size_type j, size_type size1, size_type size2) { 01805 return L::mutable_restrict1(j, i, size2, size1); 01806 } 01807 01808 static 01809 BOOST_UBLAS_INLINE 01810 size_type global_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) { 01811 return L::global_restrict2(index2, size2, index1, size1); 01812 } 01813 static 01814 BOOST_UBLAS_INLINE 01815 size_type global_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) { 01816 return L::global_restrict1(index2, size2, index1, size1); 01817 } 01818 static 01819 BOOST_UBLAS_INLINE 01820 size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) { 01821 return L::global_mutable_restrict2(index2, size2, index1, size1); 01822 } 01823 static 01824 BOOST_UBLAS_INLINE 01825 size_type global_mutable_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) { 01826 return L::global_mutable_restrict1(index2, size2, index1, size1); 01827 } 01828 }; 01829 } 01830 01831 template <class Z> 01832 struct basic_lower { 01833 typedef Z size_type; 01834 typedef lower_tag triangular_type; 01835 01836 template<class L> 01837 static 01838 BOOST_UBLAS_INLINE 01839 size_type packed_size (L, size_type size_i, size_type size_j) { 01840 return L::triangular_size (size_i, size_j); 01841 } 01842 01843 static 01844 BOOST_UBLAS_INLINE 01845 bool zero (size_type i, size_type j) { 01846 return j > i; 01847 } 01848 static 01849 BOOST_UBLAS_INLINE 01850 bool one (size_type /* i */, size_type /* j */) { 01851 return false; 01852 } 01853 static 01854 BOOST_UBLAS_INLINE 01855 bool other (size_type i, size_type j) { 01856 return j <= i; 01857 } 01858 template<class L> 01859 static 01860 BOOST_UBLAS_INLINE 01861 size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) { 01862 return L::lower_element (i, size_i, j, size_j); 01863 } 01864 01865 // return nearest valid index in column j 01866 static 01867 BOOST_UBLAS_INLINE 01868 size_type restrict1 (size_type i, size_type j, size_type size1, size_type size2) { 01869 return (std::max)(j, (std::min) (size1, i)); 01870 } 01871 // return nearest valid index in row i 01872 static 01873 BOOST_UBLAS_INLINE 01874 size_type restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) { 01875 return (std::max)(size_type(0), (std::min) (i+1, j)); 01876 } 01877 // return nearest valid mutable index in column j 01878 static 01879 BOOST_UBLAS_INLINE 01880 size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type /* size2 */) { 01881 return (std::max)(j, (std::min) (size1, i)); 01882 } 01883 // return nearest valid mutable index in row i 01884 static 01885 BOOST_UBLAS_INLINE 01886 size_type mutable_restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) { 01887 return (std::max)(size_type(0), (std::min) (i+1, j)); 01888 } 01889 01890 // return an index between the first and (1+last) filled row 01891 static 01892 BOOST_UBLAS_INLINE 01893 size_type global_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) { 01894 return (std::max)(size_type(0), (std::min)(size1, index1) ); 01895 } 01896 // return an index between the first and (1+last) filled column 01897 static 01898 BOOST_UBLAS_INLINE 01899 size_type global_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) { 01900 return (std::max)(size_type(0), (std::min)(size2, index2) ); 01901 } 01902 01903 // return an index between the first and (1+last) filled mutable row 01904 static 01905 BOOST_UBLAS_INLINE 01906 size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) { 01907 return (std::max)(size_type(0), (std::min)(size1, index1) ); 01908 } 01909 // return an index between the first and (1+last) filled mutable column 01910 static 01911 BOOST_UBLAS_INLINE 01912 size_type global_mutable_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) { 01913 return (std::max)(size_type(0), (std::min)(size2, index2) ); 01914 } 01915 }; 01916 01917 // the first row only contains a single 1. Thus it is not stored. 01918 template <class Z> 01919 struct basic_unit_lower : public basic_lower<Z> { 01920 typedef Z size_type; 01921 typedef unit_lower_tag triangular_type; 01922 01923 template<class L> 01924 static 01925 BOOST_UBLAS_INLINE 01926 size_type packed_size (L, size_type size_i, size_type size_j) { 01927 // Zero size strict triangles are bad at this point 01928 BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0, bad_index ()); 01929 return L::triangular_size (size_i - 1, size_j - 1); 01930 } 01931 01932 static 01933 BOOST_UBLAS_INLINE 01934 bool one (size_type i, size_type j) { 01935 return j == i; 01936 } 01937 static 01938 BOOST_UBLAS_INLINE 01939 bool other (size_type i, size_type j) { 01940 return j < i; 01941 } 01942 template<class L> 01943 static 01944 BOOST_UBLAS_INLINE 01945 size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) { 01946 // Zero size strict triangles are bad at this point 01947 BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0 && i != 0, bad_index ()); 01948 return L::lower_element (i-1, size_i - 1, j, size_j - 1); 01949 } 01950 01951 static 01952 BOOST_UBLAS_INLINE 01953 size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type /* size2 */) { 01954 return (std::max)(j+1, (std::min) (size1, i)); 01955 } 01956 static 01957 BOOST_UBLAS_INLINE 01958 size_type mutable_restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) { 01959 return (std::max)(size_type(0), (std::min) (i, j)); 01960 } 01961 01962 // return an index between the first and (1+last) filled mutable row 01963 static 01964 BOOST_UBLAS_INLINE 01965 size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) { 01966 return (std::max)(size_type(1), (std::min)(size1, index1) ); 01967 } 01968 // return an index between the first and (1+last) filled mutable column 01969 static 01970 BOOST_UBLAS_INLINE 01971 size_type global_mutable_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) { 01972 BOOST_UBLAS_CHECK( size2 >= 1 , external_logic() ); 01973 return (std::max)(size_type(0), (std::min)(size2-1, index2) ); 01974 } 01975 }; 01976 01977 // the first row only contains no element. Thus it is not stored. 01978 template <class Z> 01979 struct basic_strict_lower : public basic_unit_lower<Z> { 01980 typedef Z size_type; 01981 typedef strict_lower_tag triangular_type; 01982 01983 template<class L> 01984 static 01985 BOOST_UBLAS_INLINE 01986 size_type packed_size (L, size_type size_i, size_type size_j) { 01987 // Zero size strict triangles are bad at this point 01988 BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0, bad_index ()); 01989 return L::triangular_size (size_i - 1, size_j - 1); 01990 } 01991 01992 static 01993 BOOST_UBLAS_INLINE 01994 bool zero (size_type i, size_type j) { 01995 return j >= i; 01996 } 01997 static 01998 BOOST_UBLAS_INLINE 01999 bool one (size_type /*i*/, size_type /*j*/) { 02000 return false; 02001 } 02002 static 02003 BOOST_UBLAS_INLINE 02004 bool other (size_type i, size_type j) { 02005 return j < i; 02006 } 02007 template<class L> 02008 static 02009 BOOST_UBLAS_INLINE 02010 size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) { 02011 // Zero size strict triangles are bad at this point 02012 BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0 && i != 0, bad_index ()); 02013 return L::lower_element (i-1, size_i - 1, j, size_j - 1); 02014 } 02015 02016 static 02017 BOOST_UBLAS_INLINE 02018 size_type restrict1 (size_type i, size_type j, size_type size1, size_type size2) { 02019 return basic_unit_lower<Z>::mutable_restrict1(i, j, size1, size2); 02020 } 02021 static 02022 BOOST_UBLAS_INLINE 02023 size_type restrict2 (size_type i, size_type j, size_type size1, size_type size2) { 02024 return basic_unit_lower<Z>::mutable_restrict2(i, j, size1, size2); 02025 } 02026 02027 // return an index between the first and (1+last) filled row 02028 static 02029 BOOST_UBLAS_INLINE 02030 size_type global_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) { 02031 return basic_unit_lower<Z>::global_mutable_restrict1(index1, size1, index2, size2); 02032 } 02033 // return an index between the first and (1+last) filled column 02034 static 02035 BOOST_UBLAS_INLINE 02036 size_type global_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) { 02037 return basic_unit_lower<Z>::global_mutable_restrict2(index1, size1, index2, size2); 02038 } 02039 }; 02040 02041 02042 template <class Z> 02043 struct basic_upper : public detail::transposed_structure<basic_lower<Z> > 02044 { 02045 typedef upper_tag triangular_type; 02046 }; 02047 02048 template <class Z> 02049 struct basic_unit_upper : public detail::transposed_structure<basic_unit_lower<Z> > 02050 { 02051 typedef unit_upper_tag triangular_type; 02052 }; 02053 02054 template <class Z> 02055 struct basic_strict_upper : public detail::transposed_structure<basic_strict_lower<Z> > 02056 { 02057 typedef strict_upper_tag triangular_type; 02058 }; 02059 02060 02061 }}} 02062 02063 #endif