1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/os/ossrv/ossrv_pub/boost_apis/boost/numeric/ublas/functional.hpp Fri Jun 15 03:10:57 2012 +0200
1.3 @@ -0,0 +1,2087 @@
1.4 +//
1.5 +// Copyright (c) 2000-2002
1.6 +// Joerg Walter, Mathias Koch
1.7 +//
1.8 +// Permission to use, copy, modify, distribute and sell this software
1.9 +// and its documentation for any purpose is hereby granted without fee,
1.10 +// provided that the above copyright notice appear in all copies and
1.11 +// that both that copyright notice and this permission notice appear
1.12 +// in supporting documentation. The authors make no representations
1.13 +// about the suitability of this software for any purpose.
1.14 +// It is provided "as is" without express or implied warranty.
1.15 +//
1.16 +// The authors gratefully acknowledge the support of
1.17 +// GeNeSys mbH & Co. KG in producing this work.
1.18 +//
1.19 +
1.20 +#ifndef _BOOST_UBLAS_FUNCTIONAL_
1.21 +#define _BOOST_UBLAS_FUNCTIONAL_
1.22 +
1.23 +#include <functional>
1.24 +
1.25 +#include <boost/numeric/ublas/traits.hpp>
1.26 +#ifdef BOOST_UBLAS_USE_DUFF_DEVICE
1.27 +#include <boost/numeric/ublas/detail/duff.hpp>
1.28 +#endif
1.29 +#ifdef BOOST_UBLAS_USE_SIMD
1.30 +#include <boost/numeric/ublas/detail/raw.hpp>
1.31 +#else
1.32 +namespace boost { namespace numeric { namespace ublas { namespace raw {
1.33 +}}}}
1.34 +#endif
1.35 +#ifdef BOOST_UBLAS_HAVE_BINDINGS
1.36 +#include <boost/numeric/bindings/traits/std_vector.hpp>
1.37 +#include <boost/numeric/bindings/traits/ublas_vector.hpp>
1.38 +#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
1.39 +#include <boost/numeric/bindings/atlas/cblas.hpp>
1.40 +#endif
1.41 +
1.42 +#include <boost/numeric/ublas/detail/definitions.hpp>
1.43 +
1.44 +
1.45 +
1.46 +namespace boost { namespace numeric { namespace ublas {
1.47 +
1.48 + // Scalar functors
1.49 +
1.50 + // Unary
1.51 + template<class T>
1.52 + struct scalar_unary_functor {
1.53 + typedef T value_type;
1.54 + typedef typename type_traits<T>::const_reference argument_type;
1.55 + typedef typename type_traits<T>::value_type result_type;
1.56 + };
1.57 +
1.58 + template<class T>
1.59 + struct scalar_identity:
1.60 + public scalar_unary_functor<T> {
1.61 + typedef typename scalar_unary_functor<T>::argument_type argument_type;
1.62 + typedef typename scalar_unary_functor<T>::result_type result_type;
1.63 +
1.64 + static BOOST_UBLAS_INLINE
1.65 + result_type apply (argument_type t) {
1.66 + return t;
1.67 + }
1.68 + };
1.69 + template<class T>
1.70 + struct scalar_negate:
1.71 + public scalar_unary_functor<T> {
1.72 + typedef typename scalar_unary_functor<T>::argument_type argument_type;
1.73 + typedef typename scalar_unary_functor<T>::result_type result_type;
1.74 +
1.75 + static BOOST_UBLAS_INLINE
1.76 + result_type apply (argument_type t) {
1.77 + return - t;
1.78 + }
1.79 + };
1.80 + template<class T>
1.81 + struct scalar_conj:
1.82 + public scalar_unary_functor<T> {
1.83 + typedef typename scalar_unary_functor<T>::value_type value_type;
1.84 + typedef typename scalar_unary_functor<T>::argument_type argument_type;
1.85 + typedef typename scalar_unary_functor<T>::result_type result_type;
1.86 +
1.87 + static BOOST_UBLAS_INLINE
1.88 + result_type apply (argument_type t) {
1.89 + return type_traits<value_type>::conj (t);
1.90 + }
1.91 + };
1.92 +
1.93 + // Unary returning real
1.94 + template<class T>
1.95 + struct scalar_real_unary_functor {
1.96 + typedef T value_type;
1.97 + typedef typename type_traits<T>::const_reference argument_type;
1.98 + typedef typename type_traits<T>::real_type result_type;
1.99 + };
1.100 +
1.101 + template<class T>
1.102 + struct scalar_real:
1.103 + public scalar_real_unary_functor<T> {
1.104 + typedef typename scalar_real_unary_functor<T>::value_type value_type;
1.105 + typedef typename scalar_real_unary_functor<T>::argument_type argument_type;
1.106 + typedef typename scalar_real_unary_functor<T>::result_type result_type;
1.107 +
1.108 + static BOOST_UBLAS_INLINE
1.109 + result_type apply (argument_type t) {
1.110 + return type_traits<value_type>::real (t);
1.111 + }
1.112 + };
1.113 + template<class T>
1.114 + struct scalar_imag:
1.115 + public scalar_real_unary_functor<T> {
1.116 + typedef typename scalar_real_unary_functor<T>::value_type value_type;
1.117 + typedef typename scalar_real_unary_functor<T>::argument_type argument_type;
1.118 + typedef typename scalar_real_unary_functor<T>::result_type result_type;
1.119 +
1.120 + static BOOST_UBLAS_INLINE
1.121 + result_type apply (argument_type t) {
1.122 + return type_traits<value_type>::imag (t);
1.123 + }
1.124 + };
1.125 +
1.126 + // Binary
1.127 + template<class T1, class T2>
1.128 + struct scalar_binary_functor {
1.129 + typedef typename type_traits<T1>::const_reference argument1_type;
1.130 + typedef typename type_traits<T2>::const_reference argument2_type;
1.131 + typedef typename promote_traits<T1, T2>::promote_type result_type;
1.132 + };
1.133 +
1.134 + template<class T1, class T2>
1.135 + struct scalar_plus:
1.136 + public scalar_binary_functor<T1, T2> {
1.137 + typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;
1.138 + typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;
1.139 + typedef typename scalar_binary_functor<T1, T2>::result_type result_type;
1.140 +
1.141 + static BOOST_UBLAS_INLINE
1.142 + result_type apply (argument1_type t1, argument2_type t2) {
1.143 + return t1 + t2;
1.144 + }
1.145 + };
1.146 + template<class T1, class T2>
1.147 + struct scalar_minus:
1.148 + public scalar_binary_functor<T1, T2> {
1.149 + typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;
1.150 + typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;
1.151 + typedef typename scalar_binary_functor<T1, T2>::result_type result_type;
1.152 +
1.153 + static BOOST_UBLAS_INLINE
1.154 + result_type apply (argument1_type t1, argument2_type t2) {
1.155 + return t1 - t2;
1.156 + }
1.157 + };
1.158 + template<class T1, class T2>
1.159 + struct scalar_multiplies:
1.160 + public scalar_binary_functor<T1, T2> {
1.161 + typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;
1.162 + typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;
1.163 + typedef typename scalar_binary_functor<T1, T2>::result_type result_type;
1.164 +
1.165 + static BOOST_UBLAS_INLINE
1.166 + result_type apply (argument1_type t1, argument2_type t2) {
1.167 + return t1 * t2;
1.168 + }
1.169 + };
1.170 + template<class T1, class T2>
1.171 + struct scalar_divides:
1.172 + public scalar_binary_functor<T1, T2> {
1.173 + typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;
1.174 + typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;
1.175 + typedef typename scalar_binary_functor<T1, T2>::result_type result_type;
1.176 +
1.177 + static BOOST_UBLAS_INLINE
1.178 + result_type apply (argument1_type t1, argument2_type t2) {
1.179 + return t1 / t2;
1.180 + }
1.181 + };
1.182 +
1.183 + template<class T1, class T2>
1.184 + struct scalar_binary_assign_functor {
1.185 + // ISSUE Remove reference to avoid reference to reference problems
1.186 + typedef typename type_traits<typename boost::remove_reference<T1>::type>::reference argument1_type;
1.187 + typedef typename type_traits<T2>::const_reference argument2_type;
1.188 + };
1.189 +
1.190 + struct assign_tag {};
1.191 + struct computed_assign_tag {};
1.192 +
1.193 + template<class T1, class T2>
1.194 + struct scalar_assign:
1.195 + public scalar_binary_assign_functor<T1, T2> {
1.196 + typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
1.197 + typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
1.198 +#if BOOST_WORKAROUND( __IBMCPP__, <=600 )
1.199 + static const bool computed ;
1.200 +#else
1.201 + static const bool computed = false ;
1.202 +#endif
1.203 +
1.204 + static BOOST_UBLAS_INLINE
1.205 + void apply (argument1_type t1, argument2_type t2) {
1.206 + t1 = t2;
1.207 + }
1.208 +
1.209 + template<class U1, class U2>
1.210 + struct rebind {
1.211 + typedef scalar_assign<U1, U2> other;
1.212 + };
1.213 + };
1.214 +
1.215 +#if BOOST_WORKAROUND( __IBMCPP__, <=600 )
1.216 + template<class T1, class T2>
1.217 + const bool scalar_assign<T1,T2>::computed = false;
1.218 +#endif
1.219 +
1.220 + template<class T1, class T2>
1.221 + struct scalar_plus_assign:
1.222 + public scalar_binary_assign_functor<T1, T2> {
1.223 + typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
1.224 + typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
1.225 +#if BOOST_WORKAROUND( __IBMCPP__, <=600 )
1.226 + static const bool computed ;
1.227 +#else
1.228 + static const bool computed = true ;
1.229 +#endif
1.230 +
1.231 + static BOOST_UBLAS_INLINE
1.232 + void apply (argument1_type t1, argument2_type t2) {
1.233 + t1 += t2;
1.234 + }
1.235 +
1.236 + template<class U1, class U2>
1.237 + struct rebind {
1.238 + typedef scalar_plus_assign<U1, U2> other;
1.239 + };
1.240 + };
1.241 +
1.242 +#if BOOST_WORKAROUND( __IBMCPP__, <=600 )
1.243 + template<class T1, class T2>
1.244 + const bool scalar_plus_assign<T1,T2>::computed = true;
1.245 +#endif
1.246 +
1.247 + template<class T1, class T2>
1.248 + struct scalar_minus_assign:
1.249 + public scalar_binary_assign_functor<T1, T2> {
1.250 + typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
1.251 + typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
1.252 +#if BOOST_WORKAROUND( __IBMCPP__, <=600 )
1.253 + static const bool computed ;
1.254 +#else
1.255 + static const bool computed = true ;
1.256 +#endif
1.257 +
1.258 + static BOOST_UBLAS_INLINE
1.259 + void apply (argument1_type t1, argument2_type t2) {
1.260 + t1 -= t2;
1.261 + }
1.262 +
1.263 + template<class U1, class U2>
1.264 + struct rebind {
1.265 + typedef scalar_minus_assign<U1, U2> other;
1.266 + };
1.267 + };
1.268 +
1.269 +#if BOOST_WORKAROUND( __IBMCPP__, <=600 )
1.270 + template<class T1, class T2>
1.271 + const bool scalar_minus_assign<T1,T2>::computed = true;
1.272 +#endif
1.273 +
1.274 + template<class T1, class T2>
1.275 + struct scalar_multiplies_assign:
1.276 + public scalar_binary_assign_functor<T1, T2> {
1.277 + typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
1.278 + typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
1.279 + static const bool computed = true;
1.280 +
1.281 + static BOOST_UBLAS_INLINE
1.282 + void apply (argument1_type t1, argument2_type t2) {
1.283 + t1 *= t2;
1.284 + }
1.285 +
1.286 + template<class U1, class U2>
1.287 + struct rebind {
1.288 + typedef scalar_multiplies_assign<U1, U2> other;
1.289 + };
1.290 + };
1.291 + template<class T1, class T2>
1.292 + struct scalar_divides_assign:
1.293 + public scalar_binary_assign_functor<T1, T2> {
1.294 + typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
1.295 + typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
1.296 + static const bool computed ;
1.297 +
1.298 + static BOOST_UBLAS_INLINE
1.299 + void apply (argument1_type t1, argument2_type t2) {
1.300 + t1 /= t2;
1.301 + }
1.302 +
1.303 + template<class U1, class U2>
1.304 + struct rebind {
1.305 + typedef scalar_divides_assign<U1, U2> other;
1.306 + };
1.307 + };
1.308 + template<class T1, class T2>
1.309 + const bool scalar_divides_assign<T1,T2>::computed = true;
1.310 +
1.311 + template<class T1, class T2>
1.312 + struct scalar_binary_swap_functor {
1.313 + typedef typename type_traits<typename boost::remove_reference<T1>::type>::reference argument1_type;
1.314 + typedef typename type_traits<typename boost::remove_reference<T2>::type>::reference argument2_type;
1.315 + };
1.316 +
1.317 + template<class T1, class T2>
1.318 + struct scalar_swap:
1.319 + public scalar_binary_swap_functor<T1, T2> {
1.320 + typedef typename scalar_binary_swap_functor<T1, T2>::argument1_type argument1_type;
1.321 + typedef typename scalar_binary_swap_functor<T1, T2>::argument2_type argument2_type;
1.322 +
1.323 + static BOOST_UBLAS_INLINE
1.324 + void apply (argument1_type t1, argument2_type t2) {
1.325 + std::swap (t1, t2);
1.326 + }
1.327 +
1.328 + template<class U1, class U2>
1.329 + struct rebind {
1.330 + typedef scalar_swap<U1, U2> other;
1.331 + };
1.332 + };
1.333 +
1.334 + // Vector functors
1.335 +
1.336 + // Unary returning scalar
1.337 + template<class T>
1.338 + struct vector_scalar_unary_functor {
1.339 + typedef std::size_t size_type;
1.340 + typedef std::ptrdiff_t difference_type;
1.341 + typedef T value_type;
1.342 + typedef T result_type;
1.343 + };
1.344 +
1.345 + template<class T>
1.346 + struct vector_sum:
1.347 + public vector_scalar_unary_functor<T> {
1.348 + typedef typename vector_scalar_unary_functor<T>::size_type size_type;
1.349 + typedef typename vector_scalar_unary_functor<T>::difference_type difference_type;
1.350 + typedef typename vector_scalar_unary_functor<T>::value_type value_type;
1.351 + typedef typename vector_scalar_unary_functor<T>::result_type result_type;
1.352 +
1.353 + template<class E>
1.354 + static BOOST_UBLAS_INLINE
1.355 + result_type apply (const vector_expression<E> &e) {
1.356 + result_type t = result_type (0);
1.357 + size_type size (e ().size ());
1.358 + for (size_type i = 0; i < size; ++ i)
1.359 + t += e () (i);
1.360 + return t;
1.361 + }
1.362 + // Dense case
1.363 + template<class I>
1.364 + static BOOST_UBLAS_INLINE
1.365 + result_type apply (difference_type size, I it) {
1.366 + result_type t = result_type (0);
1.367 + while (-- size >= 0)
1.368 + t += *it, ++ it;
1.369 + return t;
1.370 + }
1.371 + // Sparse case
1.372 + template<class I>
1.373 + static BOOST_UBLAS_INLINE
1.374 + result_type apply (I it, const I &it_end) {
1.375 + result_type t = result_type (0);
1.376 + while (it != it_end)
1.377 + t += *it, ++ it;
1.378 + return t;
1.379 + }
1.380 + };
1.381 +
1.382 + // Unary returning real scalar
1.383 + template<class T>
1.384 + struct vector_scalar_real_unary_functor {
1.385 + typedef std::size_t size_type;
1.386 + typedef std::ptrdiff_t difference_type;
1.387 + typedef T value_type;
1.388 + typedef typename type_traits<T>::real_type real_type;
1.389 + typedef real_type result_type;
1.390 + };
1.391 +
1.392 + template<class T>
1.393 + struct vector_norm_1:
1.394 + public vector_scalar_real_unary_functor<T> {
1.395 + typedef typename vector_scalar_real_unary_functor<T>::size_type size_type;
1.396 + typedef typename vector_scalar_real_unary_functor<T>::difference_type difference_type;
1.397 + typedef typename vector_scalar_real_unary_functor<T>::value_type value_type;
1.398 + typedef typename vector_scalar_real_unary_functor<T>::real_type real_type;
1.399 + typedef typename vector_scalar_real_unary_functor<T>::result_type result_type;
1.400 +
1.401 + template<class E>
1.402 + static BOOST_UBLAS_INLINE
1.403 + result_type apply (const vector_expression<E> &e) {
1.404 + real_type t = real_type ();
1.405 + size_type size (e ().size ());
1.406 + for (size_type i = 0; i < size; ++ i) {
1.407 + real_type u (type_traits<value_type>::type_abs (e () (i)));
1.408 + t += u;
1.409 + }
1.410 + return t;
1.411 + }
1.412 + // Dense case
1.413 + template<class I>
1.414 + static BOOST_UBLAS_INLINE
1.415 + result_type apply (difference_type size, I it) {
1.416 + real_type t = real_type ();
1.417 + while (-- size >= 0) {
1.418 + real_type u (type_traits<value_type>::norm_1 (*it));
1.419 + t += u;
1.420 + ++ it;
1.421 + }
1.422 + return t;
1.423 + }
1.424 + // Sparse case
1.425 + template<class I>
1.426 + static BOOST_UBLAS_INLINE
1.427 + result_type apply (I it, const I &it_end) {
1.428 + real_type t = real_type ();
1.429 + while (it != it_end) {
1.430 + real_type u (type_traits<value_type>::norm_1 (*it));
1.431 + t += u;
1.432 + ++ it;
1.433 + }
1.434 + return t;
1.435 + }
1.436 + };
1.437 + template<class T>
1.438 + struct vector_norm_2:
1.439 + public vector_scalar_real_unary_functor<T> {
1.440 + typedef typename vector_scalar_real_unary_functor<T>::size_type size_type;
1.441 + typedef typename vector_scalar_real_unary_functor<T>::difference_type difference_type;
1.442 + typedef typename vector_scalar_real_unary_functor<T>::value_type value_type;
1.443 + typedef typename vector_scalar_real_unary_functor<T>::real_type real_type;
1.444 + typedef typename vector_scalar_real_unary_functor<T>::result_type result_type;
1.445 +
1.446 + template<class E>
1.447 + static BOOST_UBLAS_INLINE
1.448 + result_type apply (const vector_expression<E> &e) {
1.449 +#ifndef BOOST_UBLAS_SCALED_NORM
1.450 + real_type t = real_type ();
1.451 + size_type size (e ().size ());
1.452 + for (size_type i = 0; i < size; ++ i) {
1.453 + real_type u (type_traits<value_type>::norm_2 (e () (i)));
1.454 + t += u * u;
1.455 + }
1.456 + return type_traits<real_type>::type_sqrt (t);
1.457 +#else
1.458 + real_type scale = real_type ();
1.459 + real_type sum_squares (1);
1.460 + size_type size (e ().size ());
1.461 + for (size_type i = 0; i < size; ++ i) {
1.462 + real_type u (type_traits<value_type>::norm_2 (e () (i)));
1.463 + if (scale < u) {
1.464 + real_type v (scale / u);
1.465 + sum_squares = sum_squares * v * v + real_type (1);
1.466 + scale = u;
1.467 + } else {
1.468 + real_type v (u / scale);
1.469 + sum_squares += v * v;
1.470 + }
1.471 + }
1.472 + return scale * type_traits<real_type>::type_sqrt (sum_squares);
1.473 +#endif
1.474 + }
1.475 + // Dense case
1.476 + template<class I>
1.477 + static BOOST_UBLAS_INLINE
1.478 + result_type apply (difference_type size, I it) {
1.479 +#ifndef BOOST_UBLAS_SCALED_NORM
1.480 + real_type t = real_type ();
1.481 + while (-- size >= 0) {
1.482 + real_type u (type_traits<value_type>::norm_2 (*it));
1.483 + t += u * u;
1.484 + ++ it;
1.485 + }
1.486 + return type_traits<real_type>::type_sqrt (t);
1.487 +#else
1.488 + real_type scale = real_type ();
1.489 + real_type sum_squares (1);
1.490 + while (-- size >= 0) {
1.491 + real_type u (type_traits<value_type>::norm_2 (*it));
1.492 + if (scale < u) {
1.493 + real_type v (scale / u);
1.494 + sum_squares = sum_squares * v * v + real_type (1);
1.495 + scale = u;
1.496 + } else {
1.497 + real_type v (u / scale);
1.498 + sum_squares += v * v;
1.499 + }
1.500 + ++ it;
1.501 + }
1.502 + return scale * type_traits<real_type>::type_sqrt (sum_squares);
1.503 +#endif
1.504 + }
1.505 + // Sparse case
1.506 + template<class I>
1.507 + static BOOST_UBLAS_INLINE
1.508 + result_type apply (I it, const I &it_end) {
1.509 +#ifndef BOOST_UBLAS_SCALED_NORM
1.510 + real_type t = real_type ();
1.511 + while (it != it_end) {
1.512 + real_type u (type_traits<value_type>::norm_2 (*it));
1.513 + t += u * u;
1.514 + ++ it;
1.515 + }
1.516 + return type_traits<real_type>::type_sqrt (t);
1.517 +#else
1.518 + real_type scale = real_type ();
1.519 + real_type sum_squares (1);
1.520 + while (it != it_end) {
1.521 + real_type u (type_traits<value_type>::norm_2 (*it));
1.522 + if (scale < u) {
1.523 + real_type v (scale / u);
1.524 + sum_squares = sum_squares * v * v + real_type (1);
1.525 + scale = u;
1.526 + } else {
1.527 + real_type v (u / scale);
1.528 + sum_squares += v * v;
1.529 + }
1.530 + ++ it;
1.531 + }
1.532 + return scale * type_traits<real_type>::type_sqrt (sum_squares);
1.533 +#endif
1.534 + }
1.535 + };
1.536 + template<class T>
1.537 + struct vector_norm_inf:
1.538 + public vector_scalar_real_unary_functor<T> {
1.539 + typedef typename vector_scalar_real_unary_functor<T>::size_type size_type;
1.540 + typedef typename vector_scalar_real_unary_functor<T>::difference_type difference_type;
1.541 + typedef typename vector_scalar_real_unary_functor<T>::value_type value_type;
1.542 + typedef typename vector_scalar_real_unary_functor<T>::real_type real_type;
1.543 + typedef typename vector_scalar_real_unary_functor<T>::result_type result_type;
1.544 +
1.545 + template<class E>
1.546 + static BOOST_UBLAS_INLINE
1.547 + result_type apply (const vector_expression<E> &e) {
1.548 + real_type t = real_type ();
1.549 + size_type size (e ().size ());
1.550 + for (size_type i = 0; i < size; ++ i) {
1.551 + real_type u (type_traits<value_type>::norm_inf (e () (i)));
1.552 + if (u > t)
1.553 + t = u;
1.554 + }
1.555 + return t;
1.556 + }
1.557 + // Dense case
1.558 + template<class I>
1.559 + static BOOST_UBLAS_INLINE
1.560 + result_type apply (difference_type size, I it) {
1.561 + real_type t = real_type ();
1.562 + while (-- size >= 0) {
1.563 + real_type u (type_traits<value_type>::norm_inf (*it));
1.564 + if (u > t)
1.565 + t = u;
1.566 + ++ it;
1.567 + }
1.568 + return t;
1.569 + }
1.570 + // Sparse case
1.571 + template<class I>
1.572 + static BOOST_UBLAS_INLINE
1.573 + result_type apply (I it, const I &it_end) {
1.574 + real_type t = real_type ();
1.575 + while (it != it_end) {
1.576 + real_type u (type_traits<value_type>::norm_inf (*it));
1.577 + if (u > t)
1.578 + t = u;
1.579 + ++ it;
1.580 + }
1.581 + return t;
1.582 + }
1.583 + };
1.584 +
1.585 + // Unary returning index
1.586 + template<class T>
1.587 + struct vector_scalar_index_unary_functor {
1.588 + typedef std::size_t size_type;
1.589 + typedef std::ptrdiff_t difference_type;
1.590 + typedef T value_type;
1.591 + typedef typename type_traits<T>::real_type real_type;
1.592 + typedef size_type result_type;
1.593 + };
1.594 +
1.595 + template<class T>
1.596 + struct vector_index_norm_inf:
1.597 + public vector_scalar_index_unary_functor<T> {
1.598 + typedef typename vector_scalar_index_unary_functor<T>::size_type size_type;
1.599 + typedef typename vector_scalar_index_unary_functor<T>::difference_type difference_type;
1.600 + typedef typename vector_scalar_index_unary_functor<T>::value_type value_type;
1.601 + typedef typename vector_scalar_index_unary_functor<T>::real_type real_type;
1.602 + typedef typename vector_scalar_index_unary_functor<T>::result_type result_type;
1.603 +
1.604 + template<class E>
1.605 + static BOOST_UBLAS_INLINE
1.606 + result_type apply (const vector_expression<E> &e) {
1.607 + // ISSUE For CBLAS compatibility return 0 index in empty case
1.608 + result_type i_norm_inf (0);
1.609 + real_type t = real_type ();
1.610 + size_type size (e ().size ());
1.611 + for (size_type i = 0; i < size; ++ i) {
1.612 + real_type u (type_traits<value_type>::norm_inf (e () (i)));
1.613 + if (u > t) {
1.614 + i_norm_inf = i;
1.615 + t = u;
1.616 + }
1.617 + }
1.618 + return i_norm_inf;
1.619 + }
1.620 + // Dense case
1.621 + template<class I>
1.622 + static BOOST_UBLAS_INLINE
1.623 + result_type apply (difference_type size, I it) {
1.624 + // ISSUE For CBLAS compatibility return 0 index in empty case
1.625 + result_type i_norm_inf (0);
1.626 + real_type t = real_type ();
1.627 + while (-- size >= 0) {
1.628 + real_type u (type_traits<value_type>::norm_inf (*it));
1.629 + if (u > t) {
1.630 + i_norm_inf = it.index ();
1.631 + t = u;
1.632 + }
1.633 + ++ it;
1.634 + }
1.635 + return i_norm_inf;
1.636 + }
1.637 + // Sparse case
1.638 + template<class I>
1.639 + static BOOST_UBLAS_INLINE
1.640 + result_type apply (I it, const I &it_end) {
1.641 + // ISSUE For CBLAS compatibility return 0 index in empty case
1.642 + result_type i_norm_inf (0);
1.643 + real_type t = real_type ();
1.644 + while (it != it_end) {
1.645 + real_type u (type_traits<value_type>::norm_inf (*it));
1.646 + if (u > t) {
1.647 + i_norm_inf = it.index ();
1.648 + t = u;
1.649 + }
1.650 + ++ it;
1.651 + }
1.652 + return i_norm_inf;
1.653 + }
1.654 + };
1.655 +
1.656 + // Binary returning scalar
1.657 + template<class T1, class T2, class TR>
1.658 + struct vector_scalar_binary_functor {
1.659 + typedef std::size_t size_type;
1.660 + typedef std::ptrdiff_t difference_type;
1.661 + typedef TR value_type;
1.662 + typedef TR result_type;
1.663 + };
1.664 +
1.665 + template<class T1, class T2, class TR>
1.666 + struct vector_inner_prod:
1.667 + public vector_scalar_binary_functor<T1, T2, TR> {
1.668 + typedef typename vector_scalar_binary_functor<T1, T2, TR>::size_type size_type ;
1.669 + typedef typename vector_scalar_binary_functor<T1, T2, TR>::difference_type difference_type;
1.670 + typedef typename vector_scalar_binary_functor<T1, T2, TR>::value_type value_type;
1.671 + typedef typename vector_scalar_binary_functor<T1, T2, TR>::result_type result_type;
1.672 +
1.673 + template<class C1, class C2>
1.674 + static BOOST_UBLAS_INLINE
1.675 + result_type apply (const vector_container<C1> &c1,
1.676 + const vector_container<C2> &c2) {
1.677 +#ifdef BOOST_UBLAS_USE_SIMD
1.678 + using namespace raw;
1.679 + size_type size (BOOST_UBLAS_SAME (c1 ().size (), c2 ().size ()));
1.680 + const T1 *data1 = data_const (c1 ());
1.681 + const T2 *data2 = data_const (c2 ());
1.682 + size_type s1 = stride (c1 ());
1.683 + size_type s2 = stride (c2 ());
1.684 + result_type t = result_type (0);
1.685 + if (s1 == 1 && s2 == 1) {
1.686 + for (size_type i = 0; i < size; ++ i)
1.687 + t += data1 [i] * data2 [i];
1.688 + } else if (s2 == 1) {
1.689 + for (size_type i = 0, i1 = 0; i < size; ++ i, i1 += s1)
1.690 + t += data1 [i1] * data2 [i];
1.691 + } else if (s1 == 1) {
1.692 + for (size_type i = 0, i2 = 0; i < size; ++ i, i2 += s2)
1.693 + t += data1 [i] * data2 [i2];
1.694 + } else {
1.695 + for (size_type i = 0, i1 = 0, i2 = 0; i < size; ++ i, i1 += s1, i2 += s2)
1.696 + t += data1 [i1] * data2 [i2];
1.697 + }
1.698 + return t;
1.699 +#elif defined(BOOST_UBLAS_HAVE_BINDINGS)
1.700 + return boost::numeric::bindings::atlas::dot (c1 (), c2 ());
1.701 +#else
1.702 + return apply (static_cast<const vector_expression<C1> > (c1), static_cast<const vector_expression<C2> > (c2));
1.703 +#endif
1.704 + }
1.705 + template<class E1, class E2>
1.706 + static BOOST_UBLAS_INLINE
1.707 + result_type apply (const vector_expression<E1> &e1,
1.708 + const vector_expression<E2> &e2) {
1.709 + size_type size (BOOST_UBLAS_SAME (e1 ().size (), e2 ().size ()));
1.710 + result_type t = result_type (0);
1.711 +#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
1.712 + for (size_type i = 0; i < size; ++ i)
1.713 + t += e1 () (i) * e2 () (i);
1.714 +#else
1.715 + size_type i (0);
1.716 + DD (size, 4, r, (t += e1 () (i) * e2 () (i), ++ i));
1.717 +#endif
1.718 + return t;
1.719 + }
1.720 + // Dense case
1.721 + template<class I1, class I2>
1.722 + static BOOST_UBLAS_INLINE
1.723 + result_type apply (difference_type size, I1 it1, I2 it2) {
1.724 + result_type t = result_type (0);
1.725 +#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
1.726 + while (-- size >= 0)
1.727 + t += *it1 * *it2, ++ it1, ++ it2;
1.728 +#else
1.729 + DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));
1.730 +#endif
1.731 + return t;
1.732 + }
1.733 + // Packed case
1.734 + template<class I1, class I2>
1.735 + static BOOST_UBLAS_INLINE
1.736 + result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) {
1.737 + result_type t = result_type (0);
1.738 + difference_type it1_size (it1_end - it1);
1.739 + difference_type it2_size (it2_end - it2);
1.740 + difference_type diff (0);
1.741 + if (it1_size > 0 && it2_size > 0)
1.742 + diff = it2.index () - it1.index ();
1.743 + if (diff != 0) {
1.744 + difference_type size = (std::min) (diff, it1_size);
1.745 + if (size > 0) {
1.746 + it1 += size;
1.747 + it1_size -= size;
1.748 + diff -= size;
1.749 + }
1.750 + size = (std::min) (- diff, it2_size);
1.751 + if (size > 0) {
1.752 + it2 += size;
1.753 + it2_size -= size;
1.754 + diff += size;
1.755 + }
1.756 + }
1.757 + difference_type size ((std::min) (it1_size, it2_size));
1.758 + while (-- size >= 0)
1.759 + t += *it1 * *it2, ++ it1, ++ it2;
1.760 + return t;
1.761 + }
1.762 + // Sparse case
1.763 + template<class I1, class I2>
1.764 + static BOOST_UBLAS_INLINE
1.765 + result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) {
1.766 + result_type t = result_type (0);
1.767 + if (it1 != it1_end && it2 != it2_end) {
1.768 + size_type it1_index = it1.index (), it2_index = it2.index ();
1.769 + while (true) {
1.770 + difference_type compare = it1_index - it2_index;
1.771 + if (compare == 0) {
1.772 + t += *it1 * *it2, ++ it1, ++ it2;
1.773 + if (it1 != it1_end && it2 != it2_end) {
1.774 + it1_index = it1.index ();
1.775 + it2_index = it2.index ();
1.776 + } else
1.777 + break;
1.778 + } else if (compare < 0) {
1.779 + increment (it1, it1_end, - compare);
1.780 + if (it1 != it1_end)
1.781 + it1_index = it1.index ();
1.782 + else
1.783 + break;
1.784 + } else if (compare > 0) {
1.785 + increment (it2, it2_end, compare);
1.786 + if (it2 != it2_end)
1.787 + it2_index = it2.index ();
1.788 + else
1.789 + break;
1.790 + }
1.791 + }
1.792 + }
1.793 + return t;
1.794 + }
1.795 + };
1.796 +
1.797 + // Matrix functors
1.798 +
1.799 + // Binary returning vector
1.800 + template<class T1, class T2, class TR>
1.801 + struct matrix_vector_binary_functor {
1.802 + typedef std::size_t size_type;
1.803 + typedef std::ptrdiff_t difference_type;
1.804 + typedef TR value_type;
1.805 + typedef TR result_type;
1.806 + };
1.807 +
1.808 + template<class T1, class T2, class TR>
1.809 + struct matrix_vector_prod1:
1.810 + public matrix_vector_binary_functor<T1, T2, TR> {
1.811 + typedef typename matrix_vector_binary_functor<T1, T2, TR>::size_type size_type;
1.812 + typedef typename matrix_vector_binary_functor<T1, T2, TR>::difference_type difference_type;
1.813 + typedef typename matrix_vector_binary_functor<T1, T2, TR>::value_type value_type;
1.814 + typedef typename matrix_vector_binary_functor<T1, T2, TR>::result_type result_type;
1.815 +
1.816 + template<class C1, class C2>
1.817 + static BOOST_UBLAS_INLINE
1.818 + result_type apply (const matrix_container<C1> &c1,
1.819 + const vector_container<C2> &c2,
1.820 + size_type i) {
1.821 +#ifdef BOOST_UBLAS_USE_SIMD
1.822 + using namespace raw;
1.823 + size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().size ());
1.824 + const T1 *data1 = data_const (c1 ()) + i * stride1 (c1 ());
1.825 + const T2 *data2 = data_const (c2 ());
1.826 + size_type s1 = stride2 (c1 ());
1.827 + size_type s2 = stride (c2 ());
1.828 + result_type t = result_type (0);
1.829 + if (s1 == 1 && s2 == 1) {
1.830 + for (size_type j = 0; j < size; ++ j)
1.831 + t += data1 [j] * data2 [j];
1.832 + } else if (s2 == 1) {
1.833 + for (size_type j = 0, j1 = 0; j < size; ++ j, j1 += s1)
1.834 + t += data1 [j1] * data2 [j];
1.835 + } else if (s1 == 1) {
1.836 + for (size_type j = 0, j2 = 0; j < size; ++ j, j2 += s2)
1.837 + t += data1 [j] * data2 [j2];
1.838 + } else {
1.839 + for (size_type j = 0, j1 = 0, j2 = 0; j < size; ++ j, j1 += s1, j2 += s2)
1.840 + t += data1 [j1] * data2 [j2];
1.841 + }
1.842 + return t;
1.843 +#elif defined(BOOST_UBLAS_HAVE_BINDINGS)
1.844 + return boost::numeric::bindings::atlas::dot (c1 ().row (i), c2 ());
1.845 +#else
1.846 + return apply (static_cast<const matrix_expression<C1> > (c1), static_cast<const vector_expression<C2> > (c2, i));
1.847 +#endif
1.848 + }
1.849 + template<class E1, class E2>
1.850 + static BOOST_UBLAS_INLINE
1.851 + result_type apply (const matrix_expression<E1> &e1,
1.852 + const vector_expression<E2> &e2,
1.853 + size_type i) {
1.854 + size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size ());
1.855 + result_type t = result_type (0);
1.856 +#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
1.857 + for (size_type j = 0; j < size; ++ j)
1.858 + t += e1 () (i, j) * e2 () (j);
1.859 +#else
1.860 + size_type j (0);
1.861 + DD (size, 4, r, (t += e1 () (i, j) * e2 () (j), ++ j));
1.862 +#endif
1.863 + return t;
1.864 + }
1.865 + // Dense case
1.866 + template<class I1, class I2>
1.867 + static BOOST_UBLAS_INLINE
1.868 + result_type apply (difference_type size, I1 it1, I2 it2) {
1.869 + result_type t = result_type (0);
1.870 +#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
1.871 + while (-- size >= 0)
1.872 + t += *it1 * *it2, ++ it1, ++ it2;
1.873 +#else
1.874 + DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));
1.875 +#endif
1.876 + return t;
1.877 + }
1.878 + // Packed case
1.879 + template<class I1, class I2>
1.880 + static BOOST_UBLAS_INLINE
1.881 + result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) {
1.882 + result_type t = result_type (0);
1.883 + difference_type it1_size (it1_end - it1);
1.884 + difference_type it2_size (it2_end - it2);
1.885 + difference_type diff (0);
1.886 + if (it1_size > 0 && it2_size > 0)
1.887 + diff = it2.index () - it1.index2 ();
1.888 + if (diff != 0) {
1.889 + difference_type size = (std::min) (diff, it1_size);
1.890 + if (size > 0) {
1.891 + it1 += size;
1.892 + it1_size -= size;
1.893 + diff -= size;
1.894 + }
1.895 + size = (std::min) (- diff, it2_size);
1.896 + if (size > 0) {
1.897 + it2 += size;
1.898 + it2_size -= size;
1.899 + diff += size;
1.900 + }
1.901 + }
1.902 + difference_type size ((std::min) (it1_size, it2_size));
1.903 + while (-- size >= 0)
1.904 + t += *it1 * *it2, ++ it1, ++ it2;
1.905 + return t;
1.906 + }
1.907 + // Sparse case
1.908 + template<class I1, class I2>
1.909 + static BOOST_UBLAS_INLINE
1.910 + result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
1.911 + sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) {
1.912 + result_type t = result_type (0);
1.913 + if (it1 != it1_end && it2 != it2_end) {
1.914 + size_type it1_index = it1.index2 (), it2_index = it2.index ();
1.915 + while (true) {
1.916 + difference_type compare = it1_index - it2_index;
1.917 + if (compare == 0) {
1.918 + t += *it1 * *it2, ++ it1, ++ it2;
1.919 + if (it1 != it1_end && it2 != it2_end) {
1.920 + it1_index = it1.index2 ();
1.921 + it2_index = it2.index ();
1.922 + } else
1.923 + break;
1.924 + } else if (compare < 0) {
1.925 + increment (it1, it1_end, - compare);
1.926 + if (it1 != it1_end)
1.927 + it1_index = it1.index2 ();
1.928 + else
1.929 + break;
1.930 + } else if (compare > 0) {
1.931 + increment (it2, it2_end, compare);
1.932 + if (it2 != it2_end)
1.933 + it2_index = it2.index ();
1.934 + else
1.935 + break;
1.936 + }
1.937 + }
1.938 + }
1.939 + return t;
1.940 + }
1.941 + // Sparse packed case
1.942 + template<class I1, class I2>
1.943 + static BOOST_UBLAS_INLINE
1.944 + result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */,
1.945 + sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) {
1.946 + result_type t = result_type (0);
1.947 + while (it1 != it1_end) {
1.948 + t += *it1 * it2 () (it1.index2 ());
1.949 + ++ it1;
1.950 + }
1.951 + return t;
1.952 + }
1.953 + // Packed sparse case
1.954 + template<class I1, class I2>
1.955 + static BOOST_UBLAS_INLINE
1.956 + result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end,
1.957 + packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) {
1.958 + result_type t = result_type (0);
1.959 + while (it2 != it2_end) {
1.960 + t += it1 () (it1.index1 (), it2.index ()) * *it2;
1.961 + ++ it2;
1.962 + }
1.963 + return t;
1.964 + }
1.965 + // Another dispatcher
1.966 + template<class I1, class I2>
1.967 + static BOOST_UBLAS_INLINE
1.968 + result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
1.969 + sparse_bidirectional_iterator_tag) {
1.970 + typedef typename I1::iterator_category iterator1_category;
1.971 + typedef typename I2::iterator_category iterator2_category;
1.972 + return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ());
1.973 + }
1.974 + };
1.975 +
1.976 + template<class T1, class T2, class TR>
1.977 + struct matrix_vector_prod2:
1.978 + public matrix_vector_binary_functor<T1, T2, TR> {
1.979 + typedef typename matrix_vector_binary_functor<T1, T2, TR>::size_type size_type;
1.980 + typedef typename matrix_vector_binary_functor<T1, T2, TR>::difference_type difference_type;
1.981 + typedef typename matrix_vector_binary_functor<T1, T2, TR>::value_type value_type;
1.982 + typedef typename matrix_vector_binary_functor<T1, T2, TR>::result_type result_type;
1.983 +
1.984 + template<class C1, class C2>
1.985 + static BOOST_UBLAS_INLINE
1.986 + result_type apply (const vector_container<C1> &c1,
1.987 + const matrix_container<C2> &c2,
1.988 + size_type i) {
1.989 +#ifdef BOOST_UBLAS_USE_SIMD
1.990 + using namespace raw;
1.991 + size_type size = BOOST_UBLAS_SAME (c1 ().size (), c2 ().size1 ());
1.992 + const T1 *data1 = data_const (c1 ());
1.993 + const T2 *data2 = data_const (c2 ()) + i * stride2 (c2 ());
1.994 + size_type s1 = stride (c1 ());
1.995 + size_type s2 = stride1 (c2 ());
1.996 + result_type t = result_type (0);
1.997 + if (s1 == 1 && s2 == 1) {
1.998 + for (size_type j = 0; j < size; ++ j)
1.999 + t += data1 [j] * data2 [j];
1.1000 + } else if (s2 == 1) {
1.1001 + for (size_type j = 0, j1 = 0; j < size; ++ j, j1 += s1)
1.1002 + t += data1 [j1] * data2 [j];
1.1003 + } else if (s1 == 1) {
1.1004 + for (size_type j = 0, j2 = 0; j < size; ++ j, j2 += s2)
1.1005 + t += data1 [j] * data2 [j2];
1.1006 + } else {
1.1007 + for (size_type j = 0, j1 = 0, j2 = 0; j < size; ++ j, j1 += s1, j2 += s2)
1.1008 + t += data1 [j1] * data2 [j2];
1.1009 + }
1.1010 + return t;
1.1011 +#elif defined(BOOST_UBLAS_HAVE_BINDINGS)
1.1012 + return boost::numeric::bindings::atlas::dot (c1 (), c2 ().column (i));
1.1013 +#else
1.1014 + return apply (static_cast<const vector_expression<C1> > (c1), static_cast<const matrix_expression<C2> > (c2, i));
1.1015 +#endif
1.1016 + }
1.1017 + template<class E1, class E2>
1.1018 + static BOOST_UBLAS_INLINE
1.1019 + result_type apply (const vector_expression<E1> &e1,
1.1020 + const matrix_expression<E2> &e2,
1.1021 + size_type i) {
1.1022 + size_type size = BOOST_UBLAS_SAME (e1 ().size (), e2 ().size1 ());
1.1023 + result_type t = result_type (0);
1.1024 +#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
1.1025 + for (size_type j = 0; j < size; ++ j)
1.1026 + t += e1 () (j) * e2 () (j, i);
1.1027 +#else
1.1028 + size_type j (0);
1.1029 + DD (size, 4, r, (t += e1 () (j) * e2 () (j, i), ++ j));
1.1030 +#endif
1.1031 + return t;
1.1032 + }
1.1033 + // Dense case
1.1034 + template<class I1, class I2>
1.1035 + static BOOST_UBLAS_INLINE
1.1036 + result_type apply (difference_type size, I1 it1, I2 it2) {
1.1037 + result_type t = result_type (0);
1.1038 +#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
1.1039 + while (-- size >= 0)
1.1040 + t += *it1 * *it2, ++ it1, ++ it2;
1.1041 +#else
1.1042 + DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));
1.1043 +#endif
1.1044 + return t;
1.1045 + }
1.1046 + // Packed case
1.1047 + template<class I1, class I2>
1.1048 + static BOOST_UBLAS_INLINE
1.1049 + result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) {
1.1050 + result_type t = result_type (0);
1.1051 + difference_type it1_size (it1_end - it1);
1.1052 + difference_type it2_size (it2_end - it2);
1.1053 + difference_type diff (0);
1.1054 + if (it1_size > 0 && it2_size > 0)
1.1055 + diff = it2.index1 () - it1.index ();
1.1056 + if (diff != 0) {
1.1057 + difference_type size = (std::min) (diff, it1_size);
1.1058 + if (size > 0) {
1.1059 + it1 += size;
1.1060 + it1_size -= size;
1.1061 + diff -= size;
1.1062 + }
1.1063 + size = (std::min) (- diff, it2_size);
1.1064 + if (size > 0) {
1.1065 + it2 += size;
1.1066 + it2_size -= size;
1.1067 + diff += size;
1.1068 + }
1.1069 + }
1.1070 + difference_type size ((std::min) (it1_size, it2_size));
1.1071 + while (-- size >= 0)
1.1072 + t += *it1 * *it2, ++ it1, ++ it2;
1.1073 + return t;
1.1074 + }
1.1075 + // Sparse case
1.1076 + template<class I1, class I2>
1.1077 + static BOOST_UBLAS_INLINE
1.1078 + result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
1.1079 + sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) {
1.1080 + result_type t = result_type (0);
1.1081 + if (it1 != it1_end && it2 != it2_end) {
1.1082 + size_type it1_index = it1.index (), it2_index = it2.index1 ();
1.1083 + while (true) {
1.1084 + difference_type compare = it1_index - it2_index;
1.1085 + if (compare == 0) {
1.1086 + t += *it1 * *it2, ++ it1, ++ it2;
1.1087 + if (it1 != it1_end && it2 != it2_end) {
1.1088 + it1_index = it1.index ();
1.1089 + it2_index = it2.index1 ();
1.1090 + } else
1.1091 + break;
1.1092 + } else if (compare < 0) {
1.1093 + increment (it1, it1_end, - compare);
1.1094 + if (it1 != it1_end)
1.1095 + it1_index = it1.index ();
1.1096 + else
1.1097 + break;
1.1098 + } else if (compare > 0) {
1.1099 + increment (it2, it2_end, compare);
1.1100 + if (it2 != it2_end)
1.1101 + it2_index = it2.index1 ();
1.1102 + else
1.1103 + break;
1.1104 + }
1.1105 + }
1.1106 + }
1.1107 + return t;
1.1108 + }
1.1109 + // Packed sparse case
1.1110 + template<class I1, class I2>
1.1111 + static BOOST_UBLAS_INLINE
1.1112 + result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end,
1.1113 + packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) {
1.1114 + result_type t = result_type (0);
1.1115 + while (it2 != it2_end) {
1.1116 + t += it1 () (it2.index1 ()) * *it2;
1.1117 + ++ it2;
1.1118 + }
1.1119 + return t;
1.1120 + }
1.1121 + // Sparse packed case
1.1122 + template<class I1, class I2>
1.1123 + static BOOST_UBLAS_INLINE
1.1124 + result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */,
1.1125 + sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) {
1.1126 + result_type t = result_type (0);
1.1127 + while (it1 != it1_end) {
1.1128 + t += *it1 * it2 () (it1.index (), it2.index2 ());
1.1129 + ++ it1;
1.1130 + }
1.1131 + return t;
1.1132 + }
1.1133 + // Another dispatcher
1.1134 + template<class I1, class I2>
1.1135 + static BOOST_UBLAS_INLINE
1.1136 + result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
1.1137 + sparse_bidirectional_iterator_tag) {
1.1138 + typedef typename I1::iterator_category iterator1_category;
1.1139 + typedef typename I2::iterator_category iterator2_category;
1.1140 + return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ());
1.1141 + }
1.1142 + };
1.1143 +
1.1144 + // Binary returning matrix
1.1145 + template<class T1, class T2, class TR>
1.1146 + struct matrix_matrix_binary_functor {
1.1147 + typedef std::size_t size_type;
1.1148 + typedef std::ptrdiff_t difference_type;
1.1149 + typedef TR value_type;
1.1150 + typedef TR result_type;
1.1151 + };
1.1152 +
1.1153 + template<class T1, class T2, class TR>
1.1154 + struct matrix_matrix_prod:
1.1155 + public matrix_matrix_binary_functor<T1, T2, TR> {
1.1156 + typedef typename matrix_matrix_binary_functor<T1, T2, TR>::size_type size_type;
1.1157 + typedef typename matrix_matrix_binary_functor<T1, T2, TR>::difference_type difference_type;
1.1158 + typedef typename matrix_matrix_binary_functor<T1, T2, TR>::value_type value_type;
1.1159 + typedef typename matrix_matrix_binary_functor<T1, T2, TR>::result_type result_type;
1.1160 +
1.1161 + template<class C1, class C2>
1.1162 + static BOOST_UBLAS_INLINE
1.1163 + result_type apply (const matrix_container<C1> &c1,
1.1164 + const matrix_container<C2> &c2,
1.1165 + size_type i, size_type j) {
1.1166 +#ifdef BOOST_UBLAS_USE_SIMD
1.1167 + using namespace raw;
1.1168 + size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().sizc1 ());
1.1169 + const T1 *data1 = data_const (c1 ()) + i * stride1 (c1 ());
1.1170 + const T2 *data2 = data_const (c2 ()) + j * stride2 (c2 ());
1.1171 + size_type s1 = stride2 (c1 ());
1.1172 + size_type s2 = stride1 (c2 ());
1.1173 + result_type t = result_type (0);
1.1174 + if (s1 == 1 && s2 == 1) {
1.1175 + for (size_type k = 0; k < size; ++ k)
1.1176 + t += data1 [k] * data2 [k];
1.1177 + } else if (s2 == 1) {
1.1178 + for (size_type k = 0, k1 = 0; k < size; ++ k, k1 += s1)
1.1179 + t += data1 [k1] * data2 [k];
1.1180 + } else if (s1 == 1) {
1.1181 + for (size_type k = 0, k2 = 0; k < size; ++ k, k2 += s2)
1.1182 + t += data1 [k] * data2 [k2];
1.1183 + } else {
1.1184 + for (size_type k = 0, k1 = 0, k2 = 0; k < size; ++ k, k1 += s1, k2 += s2)
1.1185 + t += data1 [k1] * data2 [k2];
1.1186 + }
1.1187 + return t;
1.1188 +#elif defined(BOOST_UBLAS_HAVE_BINDINGS)
1.1189 + return boost::numeric::bindings::atlas::dot (c1 ().row (i), c2 ().column (j));
1.1190 +#else
1.1191 + return apply (static_cast<const matrix_expression<C1> > (c1), static_cast<const matrix_expression<C2> > (c2, i));
1.1192 +#endif
1.1193 + }
1.1194 + template<class E1, class E2>
1.1195 + static BOOST_UBLAS_INLINE
1.1196 + result_type apply (const matrix_expression<E1> &e1,
1.1197 + const matrix_expression<E2> &e2,
1.1198 + size_type i, size_type j) {
1.1199 + size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ());
1.1200 + result_type t = result_type (0);
1.1201 +#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
1.1202 + for (size_type k = 0; k < size; ++ k)
1.1203 + t += e1 () (i, k) * e2 () (k, j);
1.1204 +#else
1.1205 + size_type k (0);
1.1206 + DD (size, 4, r, (t += e1 () (i, k) * e2 () (k, j), ++ k));
1.1207 +#endif
1.1208 + return t;
1.1209 + }
1.1210 + // Dense case
1.1211 + template<class I1, class I2>
1.1212 + static BOOST_UBLAS_INLINE
1.1213 + result_type apply (difference_type size, I1 it1, I2 it2) {
1.1214 + result_type t = result_type (0);
1.1215 +#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
1.1216 + while (-- size >= 0)
1.1217 + t += *it1 * *it2, ++ it1, ++ it2;
1.1218 +#else
1.1219 + DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));
1.1220 +#endif
1.1221 + return t;
1.1222 + }
1.1223 + // Packed case
1.1224 + template<class I1, class I2>
1.1225 + static BOOST_UBLAS_INLINE
1.1226 + result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, packed_random_access_iterator_tag) {
1.1227 + result_type t = result_type (0);
1.1228 + difference_type it1_size (it1_end - it1);
1.1229 + difference_type it2_size (it2_end - it2);
1.1230 + difference_type diff (0);
1.1231 + if (it1_size > 0 && it2_size > 0)
1.1232 + diff = it2.index1 () - it1.index2 ();
1.1233 + if (diff != 0) {
1.1234 + difference_type size = (std::min) (diff, it1_size);
1.1235 + if (size > 0) {
1.1236 + it1 += size;
1.1237 + it1_size -= size;
1.1238 + diff -= size;
1.1239 + }
1.1240 + size = (std::min) (- diff, it2_size);
1.1241 + if (size > 0) {
1.1242 + it2 += size;
1.1243 + it2_size -= size;
1.1244 + diff += size;
1.1245 + }
1.1246 + }
1.1247 + difference_type size ((std::min) (it1_size, it2_size));
1.1248 + while (-- size >= 0)
1.1249 + t += *it1 * *it2, ++ it1, ++ it2;
1.1250 + return t;
1.1251 + }
1.1252 + // Sparse case
1.1253 + template<class I1, class I2>
1.1254 + static BOOST_UBLAS_INLINE
1.1255 + result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) {
1.1256 + result_type t = result_type (0);
1.1257 + if (it1 != it1_end && it2 != it2_end) {
1.1258 + size_type it1_index = it1.index2 (), it2_index = it2.index1 ();
1.1259 + while (true) {
1.1260 + difference_type compare = it1_index - it2_index;
1.1261 + if (compare == 0) {
1.1262 + t += *it1 * *it2, ++ it1, ++ it2;
1.1263 + if (it1 != it1_end && it2 != it2_end) {
1.1264 + it1_index = it1.index2 ();
1.1265 + it2_index = it2.index1 ();
1.1266 + } else
1.1267 + break;
1.1268 + } else if (compare < 0) {
1.1269 + increment (it1, it1_end, - compare);
1.1270 + if (it1 != it1_end)
1.1271 + it1_index = it1.index2 ();
1.1272 + else
1.1273 + break;
1.1274 + } else if (compare > 0) {
1.1275 + increment (it2, it2_end, compare);
1.1276 + if (it2 != it2_end)
1.1277 + it2_index = it2.index1 ();
1.1278 + else
1.1279 + break;
1.1280 + }
1.1281 + }
1.1282 + }
1.1283 + return t;
1.1284 + }
1.1285 + };
1.1286 +
1.1287 + // Unary returning scalar norm
1.1288 + template<class T>
1.1289 + struct matrix_scalar_real_unary_functor {
1.1290 + typedef std::size_t size_type;
1.1291 + typedef std::ptrdiff_t difference_type;
1.1292 + typedef T value_type;
1.1293 + typedef typename type_traits<T>::real_type real_type;
1.1294 + typedef real_type result_type;
1.1295 + };
1.1296 +
1.1297 + template<class T>
1.1298 + struct matrix_norm_1:
1.1299 + public matrix_scalar_real_unary_functor<T> {
1.1300 + typedef typename matrix_scalar_real_unary_functor<T>::size_type size_type;
1.1301 + typedef typename matrix_scalar_real_unary_functor<T>::difference_type difference_type;
1.1302 + typedef typename matrix_scalar_real_unary_functor<T>::value_type value_type;
1.1303 + typedef typename matrix_scalar_real_unary_functor<T>::real_type real_type;
1.1304 + typedef typename matrix_scalar_real_unary_functor<T>::result_type result_type;
1.1305 +
1.1306 + template<class E>
1.1307 + static BOOST_UBLAS_INLINE
1.1308 + result_type apply (const matrix_expression<E> &e) {
1.1309 + real_type t = real_type ();
1.1310 + size_type size2 (e ().size2 ());
1.1311 + for (size_type j = 0; j < size2; ++ j) {
1.1312 + real_type u = real_type ();
1.1313 + size_type size1 (e ().size1 ());
1.1314 + for (size_type i = 0; i < size1; ++ i) {
1.1315 + real_type v (type_traits<value_type>::norm_1 (e () (i, j)));
1.1316 + u += v;
1.1317 + }
1.1318 + if (u > t)
1.1319 + t = u;
1.1320 + }
1.1321 + return t;
1.1322 + }
1.1323 + };
1.1324 + template<class T>
1.1325 + struct matrix_norm_frobenius:
1.1326 + public matrix_scalar_real_unary_functor<T> {
1.1327 + typedef typename matrix_scalar_real_unary_functor<T>::size_type size_type;
1.1328 + typedef typename matrix_scalar_real_unary_functor<T>::difference_type difference_type;
1.1329 + typedef typename matrix_scalar_real_unary_functor<T>::value_type value_type;
1.1330 + typedef typename matrix_scalar_real_unary_functor<T>::real_type real_type;
1.1331 + typedef typename matrix_scalar_real_unary_functor<T>::result_type result_type;
1.1332 +
1.1333 + template<class E>
1.1334 + static BOOST_UBLAS_INLINE
1.1335 + result_type apply (const matrix_expression<E> &e) {
1.1336 + real_type t = real_type ();
1.1337 + size_type size1 (e ().size1 ());
1.1338 + for (size_type i = 0; i < size1; ++ i) {
1.1339 + size_type size2 (e ().size2 ());
1.1340 + for (size_type j = 0; j < size2; ++ j) {
1.1341 + real_type u (type_traits<value_type>::norm_2 (e () (i, j)));
1.1342 + t += u * u;
1.1343 + }
1.1344 + }
1.1345 + return type_traits<real_type>::type_sqrt (t);
1.1346 + }
1.1347 + };
1.1348 + template<class T>
1.1349 + struct matrix_norm_inf:
1.1350 + public matrix_scalar_real_unary_functor<T> {
1.1351 + typedef typename matrix_scalar_real_unary_functor<T>::size_type size_type;
1.1352 + typedef typename matrix_scalar_real_unary_functor<T>::difference_type difference_type;
1.1353 + typedef typename matrix_scalar_real_unary_functor<T>::value_type value_type;
1.1354 + typedef typename matrix_scalar_real_unary_functor<T>::real_type real_type;
1.1355 + typedef typename matrix_scalar_real_unary_functor<T>::result_type result_type;
1.1356 +
1.1357 + template<class E>
1.1358 + static BOOST_UBLAS_INLINE
1.1359 + result_type apply (const matrix_expression<E> &e) {
1.1360 + real_type t = real_type ();
1.1361 + size_type size1 (e ().size1 ());
1.1362 + for (size_type i = 0; i < size1; ++ i) {
1.1363 + real_type u = real_type ();
1.1364 + size_type size2 (e ().size2 ());
1.1365 + for (size_type j = 0; j < size2; ++ j) {
1.1366 + real_type v (type_traits<value_type>::norm_inf (e () (i, j)));
1.1367 + u += v;
1.1368 + }
1.1369 + if (u > t)
1.1370 + t = u;
1.1371 + }
1.1372 + return t;
1.1373 + }
1.1374 + };
1.1375 +
1.1376 + // This functor computes the address translation
1.1377 + // matrix [i] [j] -> storage [i * size2 + j]
1.1378 + template <class Z, class D>
1.1379 + struct basic_row_major {
1.1380 + typedef Z size_type;
1.1381 + typedef D difference_type;
1.1382 + typedef row_major_tag orientation_category;
1.1383 +
1.1384 + static
1.1385 + BOOST_UBLAS_INLINE
1.1386 + size_type storage_size (size_type size1, size_type size2) {
1.1387 + // Guard against size_type overflow
1.1388 + BOOST_UBLAS_CHECK (size2 == 0 || size1 <= (std::numeric_limits<size_type>::max) () / size2, bad_size ());
1.1389 + return size1 * size2;
1.1390 + }
1.1391 +
1.1392 + // Indexing
1.1393 + static
1.1394 + BOOST_UBLAS_INLINE
1.1395 + size_type element (size_type i, size_type size1, size_type j, size_type size2) {
1.1396 + BOOST_UBLAS_CHECK (i < size1, bad_index ());
1.1397 + BOOST_UBLAS_CHECK (j < size2, bad_index ());
1.1398 + detail::ignore_unused_variable_warning(size1);
1.1399 + // Guard against size_type overflow
1.1400 + BOOST_UBLAS_CHECK (i <= ((std::numeric_limits<size_type>::max) () - j) / size2, bad_index ());
1.1401 + return i * size2 + j;
1.1402 + }
1.1403 + static
1.1404 + BOOST_UBLAS_INLINE
1.1405 + size_type address (size_type i, size_type size1, size_type j, size_type size2) {
1.1406 + BOOST_UBLAS_CHECK (i <= size1, bad_index ());
1.1407 + BOOST_UBLAS_CHECK (j <= size2, bad_index ());
1.1408 + // Guard against size_type overflow - address may be size2 past end of storage
1.1409 + BOOST_UBLAS_CHECK (size2 == 0 || i <= ((std::numeric_limits<size_type>::max) () - j) / size2, bad_index ());
1.1410 + detail::ignore_unused_variable_warning(size1);
1.1411 + return i * size2 + j;
1.1412 + }
1.1413 + static
1.1414 + BOOST_UBLAS_INLINE
1.1415 + difference_type distance1 (difference_type k, size_type /* size1 */, size_type size2) {
1.1416 + return size2 != 0 ? k / size2 : 0;
1.1417 + }
1.1418 + static
1.1419 + BOOST_UBLAS_INLINE
1.1420 + difference_type distance2 (difference_type k, size_type /* size1 */, size_type /* size2 */) {
1.1421 + return k;
1.1422 + }
1.1423 + static
1.1424 + BOOST_UBLAS_INLINE
1.1425 + size_type index1 (difference_type k, size_type /* size1 */, size_type size2) {
1.1426 + return size2 != 0 ? k / size2 : 0;
1.1427 + }
1.1428 + static
1.1429 + BOOST_UBLAS_INLINE
1.1430 + size_type index2 (difference_type k, size_type /* size1 */, size_type size2) {
1.1431 + return size2 != 0 ? k % size2 : 0;
1.1432 + }
1.1433 + static
1.1434 + BOOST_UBLAS_INLINE
1.1435 + bool fast1 () {
1.1436 + return false;
1.1437 + }
1.1438 + static
1.1439 + BOOST_UBLAS_INLINE
1.1440 + size_type one1 (size_type /* size1 */, size_type size2) {
1.1441 + return size2;
1.1442 + }
1.1443 + static
1.1444 + BOOST_UBLAS_INLINE
1.1445 + bool fast2 () {
1.1446 + return true;
1.1447 + }
1.1448 + static
1.1449 + BOOST_UBLAS_INLINE
1.1450 + size_type one2 (size_type /* size1 */, size_type /* size2 */) {
1.1451 + return 1;
1.1452 + }
1.1453 +
1.1454 + static
1.1455 + BOOST_UBLAS_INLINE
1.1456 + size_type triangular_size (size_type size1, size_type size2) {
1.1457 + size_type size = (std::max) (size1, size2);
1.1458 + // Guard against size_type overflow - simplified
1.1459 + BOOST_UBLAS_CHECK (size == 0 || size / 2 < (std::numeric_limits<size_type>::max) () / size /* +1/2 */, bad_size ());
1.1460 + return ((size + 1) * size) / 2;
1.1461 + }
1.1462 + static
1.1463 + BOOST_UBLAS_INLINE
1.1464 + size_type lower_element (size_type i, size_type size1, size_type j, size_type size2) {
1.1465 + BOOST_UBLAS_CHECK (i < size1, bad_index ());
1.1466 + BOOST_UBLAS_CHECK (j < size2, bad_index ());
1.1467 + BOOST_UBLAS_CHECK (i >= j, bad_index ());
1.1468 + detail::ignore_unused_variable_warning(size1);
1.1469 + detail::ignore_unused_variable_warning(size2);
1.1470 + // FIXME size_type overflow
1.1471 + // sigma_i (i + 1) = (i + 1) * i / 2
1.1472 + // i = 0 1 2 3, sigma = 0 1 3 6
1.1473 + return ((i + 1) * i) / 2 + j;
1.1474 + }
1.1475 + static
1.1476 + BOOST_UBLAS_INLINE
1.1477 + size_type upper_element (size_type i, size_type size1, size_type j, size_type size2) {
1.1478 + BOOST_UBLAS_CHECK (i < size1, bad_index ());
1.1479 + BOOST_UBLAS_CHECK (j < size2, bad_index ());
1.1480 + BOOST_UBLAS_CHECK (i <= j, bad_index ());
1.1481 + // FIXME size_type overflow
1.1482 + // sigma_i (size - i) = size * i - i * (i - 1) / 2
1.1483 + // i = 0 1 2 3, sigma = 0 4 7 9
1.1484 + return (i * (2 * (std::max) (size1, size2) - i + 1)) / 2 + j - i;
1.1485 + }
1.1486 +
1.1487 + static
1.1488 + BOOST_UBLAS_INLINE
1.1489 + size_type element1 (size_type i, size_type size1, size_type /* j */, size_type /* size2 */) {
1.1490 + BOOST_UBLAS_CHECK (i < size1, bad_index ());
1.1491 + detail::ignore_unused_variable_warning(size1);
1.1492 + return i;
1.1493 + }
1.1494 + static
1.1495 + BOOST_UBLAS_INLINE
1.1496 + size_type element2 (size_type /* i */, size_type /* size1 */, size_type j, size_type size2) {
1.1497 + BOOST_UBLAS_CHECK (j < size2, bad_index ());
1.1498 + detail::ignore_unused_variable_warning(size2);
1.1499 + return j;
1.1500 + }
1.1501 + static
1.1502 + BOOST_UBLAS_INLINE
1.1503 + size_type address1 (size_type i, size_type size1, size_type /* j */, size_type /* size2 */) {
1.1504 + BOOST_UBLAS_CHECK (i <= size1, bad_index ());
1.1505 + detail::ignore_unused_variable_warning(size1);
1.1506 + return i;
1.1507 + }
1.1508 + static
1.1509 + BOOST_UBLAS_INLINE
1.1510 + size_type address2 (size_type /* i */, size_type /* size1 */, size_type j, size_type size2) {
1.1511 + BOOST_UBLAS_CHECK (j <= size2, bad_index ());
1.1512 + detail::ignore_unused_variable_warning(size2);
1.1513 + return j;
1.1514 + }
1.1515 + static
1.1516 + BOOST_UBLAS_INLINE
1.1517 + size_type index1 (size_type index1, size_type /* index2 */) {
1.1518 + return index1;
1.1519 + }
1.1520 + static
1.1521 + BOOST_UBLAS_INLINE
1.1522 + size_type index2 (size_type /* index1 */, size_type index2) {
1.1523 + return index2;
1.1524 + }
1.1525 + static
1.1526 + BOOST_UBLAS_INLINE
1.1527 + size_type size1 (size_type size1, size_type /* size2 */) {
1.1528 + return size1;
1.1529 + }
1.1530 + static
1.1531 + BOOST_UBLAS_INLINE
1.1532 + size_type size2 (size_type /* size1 */, size_type size2) {
1.1533 + return size2;
1.1534 + }
1.1535 +
1.1536 + // Iterating
1.1537 + template<class I>
1.1538 + static
1.1539 + BOOST_UBLAS_INLINE
1.1540 + void increment1 (I &it, size_type /* size1 */, size_type size2) {
1.1541 + it += size2;
1.1542 + }
1.1543 + template<class I>
1.1544 + static
1.1545 + BOOST_UBLAS_INLINE
1.1546 + void decrement1 (I &it, size_type /* size1 */, size_type size2) {
1.1547 + it -= size2;
1.1548 + }
1.1549 + template<class I>
1.1550 + static
1.1551 + BOOST_UBLAS_INLINE
1.1552 + void increment2 (I &it, size_type /* size1 */, size_type /* size2 */) {
1.1553 + ++ it;
1.1554 + }
1.1555 + template<class I>
1.1556 + static
1.1557 + BOOST_UBLAS_INLINE
1.1558 + void decrement2 (I &it, size_type /* size1 */, size_type /* size2 */) {
1.1559 + -- it;
1.1560 + }
1.1561 + };
1.1562 +
1.1563 + // This functor computes the address translation
1.1564 + // matrix [i] [j] -> storage [i + j * size1]
1.1565 + template <class Z, class D>
1.1566 + struct basic_column_major {
1.1567 + typedef Z size_type;
1.1568 + typedef D difference_type;
1.1569 + typedef column_major_tag orientation_category;
1.1570 +
1.1571 + static
1.1572 + BOOST_UBLAS_INLINE
1.1573 + size_type storage_size (size_type size1, size_type size2) {
1.1574 + // Guard against size_type overflow
1.1575 + BOOST_UBLAS_CHECK (size1 == 0 || size2 <= (std::numeric_limits<size_type>::max) () / size1, bad_size ());
1.1576 + return size1 * size2;
1.1577 + }
1.1578 +
1.1579 + // Indexing
1.1580 + static
1.1581 + BOOST_UBLAS_INLINE
1.1582 + size_type element (size_type i, size_type size1, size_type j, size_type size2) {
1.1583 + BOOST_UBLAS_CHECK (i < size1, bad_index ());
1.1584 + BOOST_UBLAS_CHECK (j < size2, bad_index ());
1.1585 + detail::ignore_unused_variable_warning(size2);
1.1586 + // Guard against size_type overflow
1.1587 + BOOST_UBLAS_CHECK (j <= ((std::numeric_limits<size_type>::max) () - i) / size1, bad_index ());
1.1588 + return i + j * size1;
1.1589 + }
1.1590 + static
1.1591 + BOOST_UBLAS_INLINE
1.1592 + size_type address (size_type i, size_type size1, size_type j, size_type size2) {
1.1593 + BOOST_UBLAS_CHECK (i <= size1, bad_index ());
1.1594 + BOOST_UBLAS_CHECK (j <= size2, bad_index ());
1.1595 + detail::ignore_unused_variable_warning(size2);
1.1596 + // Guard against size_type overflow - address may be size1 past end of storage
1.1597 + BOOST_UBLAS_CHECK (size1 == 0 || j <= ((std::numeric_limits<size_type>::max) () - i) / size1, bad_index ());
1.1598 + return i + j * size1;
1.1599 + }
1.1600 + static
1.1601 + BOOST_UBLAS_INLINE
1.1602 + difference_type distance1 (difference_type k, size_type /* size1 */, size_type /* size2 */) {
1.1603 + return k;
1.1604 + }
1.1605 + static
1.1606 + BOOST_UBLAS_INLINE
1.1607 + difference_type distance2 (difference_type k, size_type size1, size_type /* size2 */) {
1.1608 + return size1 != 0 ? k / size1 : 0;
1.1609 + }
1.1610 + static
1.1611 + BOOST_UBLAS_INLINE
1.1612 + size_type index1 (difference_type k, size_type size1, size_type /* size2 */) {
1.1613 + return size1 != 0 ? k % size1 : 0;
1.1614 + }
1.1615 + static
1.1616 + BOOST_UBLAS_INLINE
1.1617 + size_type index2 (difference_type k, size_type size1, size_type /* size2 */) {
1.1618 + return size1 != 0 ? k / size1 : 0;
1.1619 + }
1.1620 + static
1.1621 + BOOST_UBLAS_INLINE
1.1622 + bool fast1 () {
1.1623 + return true;
1.1624 + }
1.1625 + static
1.1626 + BOOST_UBLAS_INLINE
1.1627 + size_type one1 (size_type /* size1 */, size_type /* size2 */) {
1.1628 + return 1;
1.1629 + }
1.1630 + static
1.1631 + BOOST_UBLAS_INLINE
1.1632 + bool fast2 () {
1.1633 + return false;
1.1634 + }
1.1635 + static
1.1636 + BOOST_UBLAS_INLINE
1.1637 + size_type one2 (size_type size1, size_type /* size2 */) {
1.1638 + return size1;
1.1639 + }
1.1640 +
1.1641 + static
1.1642 + BOOST_UBLAS_INLINE
1.1643 + size_type triangular_size (size_type size1, size_type size2) {
1.1644 + size_type size = (std::max) (size1, size2);
1.1645 + // Guard against size_type overflow - simplified
1.1646 + BOOST_UBLAS_CHECK (size == 0 || size / 2 < (std::numeric_limits<size_type>::max) () / size /* +1/2 */, bad_size ());
1.1647 + return ((size + 1) * size) / 2;
1.1648 + }
1.1649 + static
1.1650 + BOOST_UBLAS_INLINE
1.1651 + size_type lower_element (size_type i, size_type size1, size_type j, size_type size2) {
1.1652 + BOOST_UBLAS_CHECK (i < size1, bad_index ());
1.1653 + BOOST_UBLAS_CHECK (j < size2, bad_index ());
1.1654 + BOOST_UBLAS_CHECK (i >= j, bad_index ());
1.1655 + // FIXME size_type overflow
1.1656 + // sigma_j (size - j) = size * j - j * (j - 1) / 2
1.1657 + // j = 0 1 2 3, sigma = 0 4 7 9
1.1658 + return i - j + (j * (2 * (std::max) (size1, size2) - j + 1)) / 2;
1.1659 + }
1.1660 + static
1.1661 + BOOST_UBLAS_INLINE
1.1662 + size_type upper_element (size_type i, size_type size1, size_type j, size_type size2) {
1.1663 + BOOST_UBLAS_CHECK (i < size1, bad_index ());
1.1664 + BOOST_UBLAS_CHECK (j < size2, bad_index ());
1.1665 + BOOST_UBLAS_CHECK (i <= j, bad_index ());
1.1666 + // FIXME size_type overflow
1.1667 + // sigma_j (j + 1) = (j + 1) * j / 2
1.1668 + // j = 0 1 2 3, sigma = 0 1 3 6
1.1669 + return i + ((j + 1) * j) / 2;
1.1670 + }
1.1671 +
1.1672 + static
1.1673 + BOOST_UBLAS_INLINE
1.1674 + size_type element1 (size_type /* i */, size_type /* size1 */, size_type j, size_type size2) {
1.1675 + BOOST_UBLAS_CHECK (j < size2, bad_index ());
1.1676 + detail::ignore_unused_variable_warning(size2);
1.1677 + return j;
1.1678 + }
1.1679 + static
1.1680 + BOOST_UBLAS_INLINE
1.1681 + size_type element2 (size_type i, size_type size1, size_type /* j */, size_type /* size2 */) {
1.1682 + BOOST_UBLAS_CHECK (i < size1, bad_index ());
1.1683 + detail::ignore_unused_variable_warning(size1);
1.1684 + return i;
1.1685 + }
1.1686 + static
1.1687 + BOOST_UBLAS_INLINE
1.1688 + size_type address1 (size_type /* i */, size_type /* size1 */, size_type j, size_type size2) {
1.1689 + BOOST_UBLAS_CHECK (j <= size2, bad_index ());
1.1690 + detail::ignore_unused_variable_warning(size2);
1.1691 + return j;
1.1692 + }
1.1693 + static
1.1694 + BOOST_UBLAS_INLINE
1.1695 + size_type address2 (size_type i, size_type size1, size_type /* j */, size_type /* size2 */) {
1.1696 + BOOST_UBLAS_CHECK (i <= size1, bad_index ());
1.1697 + detail::ignore_unused_variable_warning(size1);
1.1698 + return i;
1.1699 + }
1.1700 + static
1.1701 + BOOST_UBLAS_INLINE
1.1702 + size_type index1 (size_type /* index1 */, size_type index2) {
1.1703 + return index2;
1.1704 + }
1.1705 + static
1.1706 + BOOST_UBLAS_INLINE
1.1707 + size_type index2 (size_type index1, size_type /* index2 */) {
1.1708 + return index1;
1.1709 + }
1.1710 + static
1.1711 + BOOST_UBLAS_INLINE
1.1712 + size_type size1 (size_type /* size1 */, size_type size2) {
1.1713 + return size2;
1.1714 + }
1.1715 + static
1.1716 + BOOST_UBLAS_INLINE
1.1717 + size_type size2 (size_type size1, size_type /* size2 */) {
1.1718 + return size1;
1.1719 + }
1.1720 +
1.1721 + // Iterating
1.1722 + template<class I>
1.1723 + static
1.1724 + BOOST_UBLAS_INLINE
1.1725 + void increment1 (I &it, size_type /* size1 */, size_type /* size2 */) {
1.1726 + ++ it;
1.1727 + }
1.1728 + template<class I>
1.1729 + static
1.1730 + BOOST_UBLAS_INLINE
1.1731 + void decrement1 (I &it, size_type /* size1 */, size_type /* size2 */) {
1.1732 + -- it;
1.1733 + }
1.1734 + template<class I>
1.1735 + static
1.1736 + BOOST_UBLAS_INLINE
1.1737 + void increment2 (I &it, size_type size1, size_type /* size2 */) {
1.1738 + it += size1;
1.1739 + }
1.1740 + template<class I>
1.1741 + static
1.1742 + BOOST_UBLAS_INLINE
1.1743 + void decrement2 (I &it, size_type size1, size_type /* size2 */) {
1.1744 + it -= size1;
1.1745 + }
1.1746 + };
1.1747 +
1.1748 +
1.1749 + template <class Z>
1.1750 + struct basic_full {
1.1751 + typedef Z size_type;
1.1752 +
1.1753 + template<class L>
1.1754 + static
1.1755 + BOOST_UBLAS_INLINE
1.1756 + size_type packed_size (size_type size1, size_type size2) {
1.1757 + return L::storage_size (size1, size2);
1.1758 + }
1.1759 +
1.1760 + static
1.1761 + BOOST_UBLAS_INLINE
1.1762 + bool zero (size_type /* i */, size_type /* j */) {
1.1763 + return false;
1.1764 + }
1.1765 + static
1.1766 + BOOST_UBLAS_INLINE
1.1767 + bool one (size_type /* i */, size_type /* j */) {
1.1768 + return false;
1.1769 + }
1.1770 + static
1.1771 + BOOST_UBLAS_INLINE
1.1772 + bool other (size_type /* i */, size_type /* j */) {
1.1773 + return true;
1.1774 + }
1.1775 + };
1.1776 +
1.1777 + template <class Z>
1.1778 + struct basic_lower {
1.1779 + typedef Z size_type;
1.1780 +
1.1781 + template<class L>
1.1782 + static
1.1783 + BOOST_UBLAS_INLINE
1.1784 + size_type packed_size (L, size_type size1, size_type size2) {
1.1785 + return L::triangular_size (size1, size2);
1.1786 + }
1.1787 +
1.1788 + static
1.1789 + BOOST_UBLAS_INLINE
1.1790 + bool zero (size_type i, size_type j) {
1.1791 + return j > i;
1.1792 + }
1.1793 + static
1.1794 + BOOST_UBLAS_INLINE
1.1795 + bool one (size_type /* i */, size_type /* j */) {
1.1796 + return false;
1.1797 + }
1.1798 + static
1.1799 + BOOST_UBLAS_INLINE
1.1800 + bool other (size_type i, size_type j) {
1.1801 + return j <= i;
1.1802 + }
1.1803 + template<class L>
1.1804 + static
1.1805 + BOOST_UBLAS_INLINE
1.1806 + size_type element (L, size_type i, size_type size1, size_type j, size_type size2) {
1.1807 + return L::lower_element (i, size1, j, size2);
1.1808 + }
1.1809 +
1.1810 + static
1.1811 + BOOST_UBLAS_INLINE
1.1812 + size_type restrict1 (size_type i, size_type j) {
1.1813 + return (std::max) (i, j);
1.1814 + }
1.1815 + static
1.1816 + BOOST_UBLAS_INLINE
1.1817 + size_type restrict2 (size_type i, size_type j) {
1.1818 + return (std::min) (i + 1, j);
1.1819 + }
1.1820 + static
1.1821 + BOOST_UBLAS_INLINE
1.1822 + size_type mutable_restrict1 (size_type i, size_type j) {
1.1823 + return (std::max) (i, j);
1.1824 + }
1.1825 + static
1.1826 + BOOST_UBLAS_INLINE
1.1827 + size_type mutable_restrict2 (size_type i, size_type j) {
1.1828 + return (std::min) (i + 1, j);
1.1829 + }
1.1830 + };
1.1831 + template <class Z>
1.1832 + struct basic_upper {
1.1833 + typedef Z size_type;
1.1834 +
1.1835 + template<class L>
1.1836 + static
1.1837 + BOOST_UBLAS_INLINE
1.1838 + size_type packed_size (L, size_type size1, size_type size2) {
1.1839 + return L::triangular_size (size1, size2);
1.1840 + }
1.1841 +
1.1842 + static
1.1843 + BOOST_UBLAS_INLINE
1.1844 + bool zero (size_type i, size_type j) {
1.1845 + return j < i;
1.1846 + }
1.1847 + static
1.1848 + BOOST_UBLAS_INLINE
1.1849 + bool one (size_type /* i */, size_type /* j */) {
1.1850 + return false;
1.1851 + }
1.1852 + static
1.1853 + BOOST_UBLAS_INLINE
1.1854 + bool other (size_type i, size_type j) {
1.1855 + return j >= i;
1.1856 + }
1.1857 + template<class L>
1.1858 + static
1.1859 + BOOST_UBLAS_INLINE
1.1860 + size_type element (L, size_type i, size_type size1, size_type j, size_type size2) {
1.1861 + return L::upper_element (i, size1, j, size2);
1.1862 + }
1.1863 +
1.1864 + static
1.1865 + BOOST_UBLAS_INLINE
1.1866 + size_type restrict1 (size_type i, size_type j) {
1.1867 + return (std::min) (i, j + 1);
1.1868 + }
1.1869 + static
1.1870 + BOOST_UBLAS_INLINE
1.1871 + size_type restrict2 (size_type i, size_type j) {
1.1872 + return (std::max) (i, j);
1.1873 + }
1.1874 + static
1.1875 + BOOST_UBLAS_INLINE
1.1876 + size_type mutable_restrict1 (size_type i, size_type j) {
1.1877 + return (std::min) (i, j + 1);
1.1878 + }
1.1879 + static
1.1880 + BOOST_UBLAS_INLINE
1.1881 + size_type mutable_restrict2 (size_type i, size_type j) {
1.1882 + return (std::max) (i, j);
1.1883 + }
1.1884 + };
1.1885 + template <class Z>
1.1886 + struct basic_unit_lower : public basic_lower<Z> {
1.1887 + typedef Z size_type;
1.1888 +
1.1889 + template<class L>
1.1890 + static
1.1891 + BOOST_UBLAS_INLINE
1.1892 + size_type packed_size (L, size_type size1, size_type size2) {
1.1893 + // Zero size strict triangles are bad at this point
1.1894 + BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());
1.1895 + return L::triangular_size (size1 - 1, size2 - 1);
1.1896 + }
1.1897 +
1.1898 + static
1.1899 + BOOST_UBLAS_INLINE
1.1900 + bool one (size_type i, size_type j) {
1.1901 + return j == i;
1.1902 + }
1.1903 + static
1.1904 + BOOST_UBLAS_INLINE
1.1905 + bool other (size_type i, size_type j) {
1.1906 + return j < i;
1.1907 + }
1.1908 + template<class L>
1.1909 + static
1.1910 + BOOST_UBLAS_INLINE
1.1911 + size_type element (L, size_type i, size_type size1, size_type j, size_type size2) {
1.1912 + // Zero size strict triangles are bad at this point
1.1913 + BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());
1.1914 + return L::lower_element (i, size1 - 1, j, size2 - 1);
1.1915 + }
1.1916 +
1.1917 + static
1.1918 + BOOST_UBLAS_INLINE
1.1919 + size_type mutable_restrict1 (size_type i, size_type j) {
1.1920 + return (std::max) (i, j);
1.1921 + }
1.1922 + static
1.1923 + BOOST_UBLAS_INLINE
1.1924 + size_type mutable_restrict2 (size_type i, size_type j) {
1.1925 + return (std::min) (i, j);
1.1926 + }
1.1927 + };
1.1928 + template <class Z>
1.1929 + struct basic_unit_upper : public basic_upper<Z> {
1.1930 + typedef Z size_type;
1.1931 +
1.1932 + template<class L>
1.1933 + static
1.1934 + BOOST_UBLAS_INLINE
1.1935 + size_type packed_size (L, size_type size1, size_type size2) {
1.1936 + // Zero size strict triangles are bad at this point
1.1937 + BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());
1.1938 + return L::triangular_size (size1 - 1, size2 - 1);
1.1939 + }
1.1940 +
1.1941 + static
1.1942 + BOOST_UBLAS_INLINE
1.1943 + bool one (size_type i, size_type j) {
1.1944 + return j == i;
1.1945 + }
1.1946 + static
1.1947 + BOOST_UBLAS_INLINE
1.1948 + bool other (size_type i, size_type j) {
1.1949 + return j > i;
1.1950 + }
1.1951 + template<class L>
1.1952 + static
1.1953 + BOOST_UBLAS_INLINE
1.1954 + size_type element (L, size_type i, size_type size1, size_type j, size_type size2) {
1.1955 + // Zero size strict triangles are bad at this point
1.1956 + BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());
1.1957 + return L::upper_element (i, size1 - 1, j, size2 - 1);
1.1958 + }
1.1959 +
1.1960 + static
1.1961 + BOOST_UBLAS_INLINE
1.1962 + size_type mutable_restrict1 (size_type i, size_type j) {
1.1963 + return (std::min) (i, j);
1.1964 + }
1.1965 + static
1.1966 + BOOST_UBLAS_INLINE
1.1967 + size_type mutable_restrict2 (size_type i, size_type j) {
1.1968 + return (std::max) (i, j);
1.1969 + }
1.1970 + };
1.1971 + template <class Z>
1.1972 + struct basic_strict_lower : public basic_lower<Z> {
1.1973 + typedef Z size_type;
1.1974 +
1.1975 + template<class L>
1.1976 + static
1.1977 + BOOST_UBLAS_INLINE
1.1978 + size_type packed_size (L, size_type size1, size_type size2) {
1.1979 + // Zero size strict triangles are bad at this point
1.1980 + BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());
1.1981 + return L::triangular_size (size1 - 1, size2 - 1);
1.1982 + }
1.1983 +
1.1984 + static
1.1985 + BOOST_UBLAS_INLINE
1.1986 + bool zero (size_type i, size_type j) {
1.1987 + return j >= i;
1.1988 + }
1.1989 + static
1.1990 + BOOST_UBLAS_INLINE
1.1991 + bool one (size_type /*i*/, size_type /*j*/) {
1.1992 + return false;
1.1993 + }
1.1994 + static
1.1995 + BOOST_UBLAS_INLINE
1.1996 + bool other (size_type i, size_type j) {
1.1997 + return j < i;
1.1998 + }
1.1999 + template<class L>
1.2000 + static
1.2001 + BOOST_UBLAS_INLINE
1.2002 + size_type element (L, size_type i, size_type size1, size_type j, size_type size2) {
1.2003 + // Zero size strict triangles are bad at this point
1.2004 + BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());
1.2005 + return L::lower_element (i, size1 - 1, j, size2 - 1);
1.2006 + }
1.2007 +
1.2008 + static
1.2009 + BOOST_UBLAS_INLINE
1.2010 + size_type restrict1 (size_type i, size_type j) {
1.2011 + return (std::max) (i, j);
1.2012 + }
1.2013 + static
1.2014 + BOOST_UBLAS_INLINE
1.2015 + size_type restrict2 (size_type i, size_type j) {
1.2016 + return (std::min) (i, j);
1.2017 + }
1.2018 + static
1.2019 + BOOST_UBLAS_INLINE
1.2020 + size_type mutable_restrict1 (size_type i, size_type j) {
1.2021 + return (std::max) (i, j);
1.2022 + }
1.2023 + static
1.2024 + BOOST_UBLAS_INLINE
1.2025 + size_type mutable_restrict2 (size_type i, size_type j) {
1.2026 + return (std::min) (i, j);
1.2027 + }
1.2028 + };
1.2029 + template <class Z>
1.2030 + struct basic_strict_upper : public basic_upper<Z> {
1.2031 + typedef Z size_type;
1.2032 +
1.2033 + template<class L>
1.2034 + static
1.2035 + BOOST_UBLAS_INLINE
1.2036 + size_type packed_size (L, size_type size1, size_type size2) {
1.2037 + // Zero size strict triangles are bad at this point
1.2038 + BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());
1.2039 + return L::triangular_size (size1 - 1, size2 - 1);
1.2040 + }
1.2041 +
1.2042 + static
1.2043 + BOOST_UBLAS_INLINE
1.2044 + bool zero (size_type i, size_type j) {
1.2045 + return j <= i;
1.2046 + }
1.2047 + static
1.2048 + BOOST_UBLAS_INLINE
1.2049 + bool one (size_type /*i*/, size_type /*j*/) {
1.2050 + return false;
1.2051 + }
1.2052 + static
1.2053 + BOOST_UBLAS_INLINE
1.2054 + bool other (size_type i, size_type j) {
1.2055 + return j > i;
1.2056 + }
1.2057 + template<class L>
1.2058 + static
1.2059 + BOOST_UBLAS_INLINE
1.2060 + size_type element (L, size_type i, size_type size1, size_type j, size_type size2) {
1.2061 + // Zero size strict triangles are bad at this point
1.2062 + BOOST_UBLAS_CHECK (size1 != 0 && size2 != 0, bad_index ());
1.2063 + return L::upper_element (i, size1 - 1, j, size2 - 1);
1.2064 + }
1.2065 +
1.2066 + static
1.2067 + BOOST_UBLAS_INLINE
1.2068 + size_type restrict1 (size_type i, size_type j) {
1.2069 + return (std::min) (i, j);
1.2070 + }
1.2071 + static
1.2072 + BOOST_UBLAS_INLINE
1.2073 + size_type restrict2 (size_type i, size_type j) {
1.2074 + return (std::max) (i, j);
1.2075 + }
1.2076 + static
1.2077 + BOOST_UBLAS_INLINE
1.2078 + size_type mutable_restrict1 (size_type i, size_type j) {
1.2079 + return (std::min) (i, j);
1.2080 + }
1.2081 + static
1.2082 + BOOST_UBLAS_INLINE
1.2083 + size_type mutable_restrict2 (size_type i, size_type j) {
1.2084 + return (std::max) (i, j);
1.2085 + }
1.2086 + };
1.2087 +
1.2088 +}}}
1.2089 +
1.2090 +#endif