00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
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
00042
00043
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
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
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
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
00328
00329
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
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
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
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
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
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 () == 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
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
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
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
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
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
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
00604 template<class D, class I>
00605 static BOOST_UBLAS_INLINE
00606 result_type apply (D size, I it) {
00607
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
00621 template<class I>
00622 static BOOST_UBLAS_INLINE
00623 result_type apply (I it, const I &it_end) {
00624
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
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
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
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
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
00771
00772
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
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
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
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
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 &,
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
00927 template<class I1, class I2>
00928 static BOOST_UBLAS_INLINE
00929 result_type apply (I1 it1, const I1 &, 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
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
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
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
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
01083 template<class I1, class I2>
01084 static BOOST_UBLAS_INLINE
01085 result_type apply (I1 it1, const I1 &, 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
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 &,
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
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
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
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
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
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
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
01347 template <class Z, class D> struct basic_column_major;
01348
01349
01350
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
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
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
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
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
01389 static
01390 BOOST_UBLAS_INLINE
01391 difference_type distance_i (difference_type k, size_type , 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_type ) {
01397 return k;
01398 }
01399 static
01400 BOOST_UBLAS_INLINE
01401 size_type index_i (difference_type k, size_type , 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_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
01421 template<class I>
01422 static
01423 BOOST_UBLAS_INLINE
01424 void increment_i (I &it, size_type , 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_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_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_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_type ) {
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_type ) {
01455 it += n;
01456 }
01457 template<class I>
01458 static
01459 BOOST_UBLAS_INLINE
01460 void decrement_j (I &it, size_type , size_type ) {
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_type ) {
01467 it -= n;
01468 }
01469
01470
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
01476 BOOST_UBLAS_CHECK (size == 0 || size / 2 < (std::numeric_limits<size_type>::max) () / size , 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
01488
01489
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
01499
01500
01501 return (i * (2 * (std::max) (size_i, size_j) - i + 1)) / 2 + j - i;
01502 }
01503
01504
01505 static
01506 BOOST_UBLAS_INLINE
01507 size_type index_M (size_type index1, size_type ) {
01508 return index1;
01509 }
01510 static
01511 BOOST_UBLAS_INLINE
01512 size_type index_m (size_type , size_type index2) {
01513 return index2;
01514 }
01515 static
01516 BOOST_UBLAS_INLINE
01517 size_type size_M (size_type size_i, size_type ) {
01518 return size_i;
01519 }
01520 static
01521 BOOST_UBLAS_INLINE
01522 size_type size_m (size_type , size_type size_j) {
01523 return size_j;
01524 }
01525 };
01526
01527
01528
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
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
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
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
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
01567 static
01568 BOOST_UBLAS_INLINE
01569 difference_type distance_i (difference_type k, size_type , size_type ) {
01570 return k;
01571 }
01572 static
01573 BOOST_UBLAS_INLINE
01574 difference_type distance_j (difference_type k, size_type size_i, size_type ) {
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 ) {
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 ) {
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
01599 template<class I>
01600 static
01601 BOOST_UBLAS_INLINE
01602 void increment_i (I &it, size_type , size_type ) {
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_type ) {
01609 it += n;
01610 }
01611 template<class I>
01612 static
01613 BOOST_UBLAS_INLINE
01614 void decrement_i (I &it, size_type , size_type ) {
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_type ) {
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 ) {
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 ) {
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 ) {
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 ) {
01645 it -= n* size_i;
01646 }
01647
01648
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
01654 BOOST_UBLAS_CHECK (size == 0 || size / 2 < (std::numeric_limits<size_type>::max) () / size , 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
01664
01665
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
01675
01676
01677 return i + ((j + 1) * j) / 2;
01678 }
01679
01680
01681 static
01682 BOOST_UBLAS_INLINE
01683 size_type index_M (size_type , size_type index2) {
01684 return index2;
01685 }
01686 static
01687 BOOST_UBLAS_INLINE
01688 size_type index_m (size_type index1, size_type ) {
01689 return index1;
01690 }
01691 static
01692 BOOST_UBLAS_INLINE
01693 size_type size_M (size_type , 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 ) {
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 , size_type ) {
01718 return false;
01719 }
01720 static
01721 BOOST_UBLAS_INLINE
01722 bool one (size_type , size_type ) {
01723 return false;
01724 }
01725 static
01726 BOOST_UBLAS_INLINE
01727 bool other (size_type , size_type ) {
01728 return true;
01729 }
01730
01731 static
01732 BOOST_UBLAS_INLINE
01733 size_type restrict1 (size_type i, size_type ) {
01734 return i;
01735 }
01736 static
01737 BOOST_UBLAS_INLINE
01738 size_type restrict2 (size_type , size_type j) {
01739 return j;
01740 }
01741 static
01742 BOOST_UBLAS_INLINE
01743 size_type mutable_restrict1 (size_type i, size_type ) {
01744 return i;
01745 }
01746 static
01747 BOOST_UBLAS_INLINE
01748 size_type mutable_restrict2 (size_type , 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 , 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 , size_type ) {
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
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
01872 static
01873 BOOST_UBLAS_INLINE
01874 size_type restrict2 (size_type i, size_type j, size_type , size_type ) {
01875 return (std::max)(size_type(0), (std::min) (i+1, j));
01876 }
01877
01878 static
01879 BOOST_UBLAS_INLINE
01880 size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type ) {
01881 return (std::max)(j, (std::min) (size1, i));
01882 }
01883
01884 static
01885 BOOST_UBLAS_INLINE
01886 size_type mutable_restrict2 (size_type i, size_type j, size_type , size_type ) {
01887 return (std::max)(size_type(0), (std::min) (i+1, j));
01888 }
01889
01890
01891 static
01892 BOOST_UBLAS_INLINE
01893 size_type global_restrict1 (size_type index1, size_type size1, size_type , size_type ) {
01894 return (std::max)(size_type(0), (std::min)(size1, index1) );
01895 }
01896
01897 static
01898 BOOST_UBLAS_INLINE
01899 size_type global_restrict2 (size_type , size_type , size_type index2, size_type size2) {
01900 return (std::max)(size_type(0), (std::min)(size2, index2) );
01901 }
01902
01903
01904 static
01905 BOOST_UBLAS_INLINE
01906 size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type , size_type ) {
01907 return (std::max)(size_type(0), (std::min)(size1, index1) );
01908 }
01909
01910 static
01911 BOOST_UBLAS_INLINE
01912 size_type global_mutable_restrict2 (size_type , size_type , size_type index2, size_type size2) {
01913 return (std::max)(size_type(0), (std::min)(size2, index2) );
01914 }
01915 };
01916
01917
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
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
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 ) {
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 , size_type ) {
01959 return (std::max)(size_type(0), (std::min) (i, j));
01960 }
01961
01962
01963 static
01964 BOOST_UBLAS_INLINE
01965 size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type , size_type ) {
01966 return (std::max)(size_type(1), (std::min)(size1, index1) );
01967 }
01968
01969 static
01970 BOOST_UBLAS_INLINE
01971 size_type global_mutable_restrict2 (size_type , size_type , 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
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
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 , size_type ) {
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
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
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
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