ceph/src/boost/libs/numeric/ublas/include/boost/numeric/ublas/functional.hpp

   1 //
   2 //  Copyright (c) 2000-2009
   3 //  Joerg Walter, Mathias Koch, Gunter Winkler
   4 //
   5 //  Distributed under the Boost Software License, Version 1.0. (See
   6 //  accompanying file LICENSE_1_0.txt or copy at
   7 //  http://www.boost.org/LICENSE_1_0.txt)
   8 //
   9 //  The authors gratefully acknowledge the support of
  10 //  GeNeSys mbH & Co. KG in producing this work.
  11 //
  12
  13 #ifndef _BOOST_UBLAS_FUNCTIONAL_
  14 #define _BOOST_UBLAS_FUNCTIONAL_
  15
  16 #include <functional>
  17
  18 #include <boost/core/ignore_unused.hpp>
  19
  20 #include <boost/numeric/ublas/traits.hpp>
  21 #ifdef BOOST_UBLAS_USE_DUFF_DEVICE
  22 #include <boost/numeric/ublas/detail/duff.hpp>
  23 #endif
  24 #ifdef BOOST_UBLAS_USE_SIMD
  25 #include <boost/numeric/ublas/detail/raw.hpp>
  26 #else
  27 namespace boost { namespace numeric { namespace ublas { namespace raw {
  28 }}}}
  29 #endif
  30 #ifdef BOOST_UBLAS_HAVE_BINDINGS
  31 #include <boost/numeric/bindings/traits/std_vector.hpp>
  32 #include <boost/numeric/bindings/traits/ublas_vector.hpp>
  33 #include <boost/numeric/bindings/traits/ublas_matrix.hpp>
  34 #include <boost/numeric/bindings/atlas/cblas.hpp>
  35 #endif
  36
  37 #include <boost/numeric/ublas/detail/definitions.hpp>
  38
  39
  40
  41 namespace boost { namespace numeric { namespace ublas {
  42
  43     // Scalar functors
  44
  45     // Unary
  46     template<class T>
  47     struct scalar_unary_functor {
  48         typedef T value_type;
  49         typedef typename type_traits<T>::const_reference argument_type;
  50         typedef typename type_traits<T>::value_type result_type;
  51     };
  52
  53     template<class T>
  54     struct scalar_identity:
  55         public scalar_unary_functor<T> {
  56         typedef typename scalar_unary_functor<T>::argument_type argument_type;
  57         typedef typename scalar_unary_functor<T>::result_type result_type;
  58
  59         static BOOST_UBLAS_INLINE
  60         result_type apply (argument_type t) {
  61             return t;
  62         }
  63     };
  64     template<class T>
  65     struct scalar_negate:
  66         public scalar_unary_functor<T> {
  67         typedef typename scalar_unary_functor<T>::argument_type argument_type;
  68         typedef typename scalar_unary_functor<T>::result_type result_type;
  69
  70         static BOOST_UBLAS_INLINE
  71         result_type apply (argument_type t) {
  72             return - t;
  73         }
  74     };
  75     template<class T>
  76     struct scalar_conj:
  77         public scalar_unary_functor<T> {
  78         typedef typename scalar_unary_functor<T>::value_type value_type;
  79         typedef typename scalar_unary_functor<T>::argument_type argument_type;
  80         typedef typename scalar_unary_functor<T>::result_type result_type;
  81
  82         static BOOST_UBLAS_INLINE
  83         result_type apply (argument_type t) {
  84             return type_traits<value_type>::conj (t);
  85         }
  86     };
  87
  88     // Unary returning real
  89     template<class T>
  90     struct scalar_real_unary_functor {
  91         typedef T value_type;
  92         typedef typename type_traits<T>::const_reference argument_type;
  93         typedef typename type_traits<T>::real_type result_type;
  94     };
  95
  96     template<class T>
  97     struct scalar_real:
  98         public scalar_real_unary_functor<T> {
  99         typedef typename scalar_real_unary_functor<T>::value_type value_type;
 100         typedef typename scalar_real_unary_functor<T>::argument_type argument_type;
 101         typedef typename scalar_real_unary_functor<T>::result_type result_type;
 102
 103         static BOOST_UBLAS_INLINE
 104         result_type apply (argument_type t) {
 105             return type_traits<value_type>::real (t);
 106         }
 107     };
 108     template<class T>
 109     struct scalar_imag:
 110         public scalar_real_unary_functor<T> {
 111         typedef typename scalar_real_unary_functor<T>::value_type value_type;
 112         typedef typename scalar_real_unary_functor<T>::argument_type argument_type;
 113         typedef typename scalar_real_unary_functor<T>::result_type result_type;
 114
 115         static BOOST_UBLAS_INLINE
 116         result_type apply (argument_type t) {
 117             return type_traits<value_type>::imag (t);
 118         }
 119     };
 120
 121     // Binary
 122     template<class T1, class T2>
 123     struct scalar_binary_functor {
 124         typedef typename type_traits<T1>::const_reference argument1_type;
 125         typedef typename type_traits<T2>::const_reference argument2_type;
 126         typedef typename promote_traits<T1, T2>::promote_type result_type;
 127     };
 128
 129     template<class T1, class T2>
 130     struct scalar_plus:
 131         public scalar_binary_functor<T1, T2> {
 132         typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;
 133         typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;
 134         typedef typename scalar_binary_functor<T1, T2>::result_type result_type;
 135
 136         static BOOST_UBLAS_INLINE
 137         result_type apply (argument1_type t1, argument2_type t2) {
 138             return t1 + t2;
 139         }
 140     };
 141     template<class T1, class T2>
 142     struct scalar_minus:
 143         public scalar_binary_functor<T1, T2> {
 144         typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;
 145         typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;
 146         typedef typename scalar_binary_functor<T1, T2>::result_type result_type;
 147
 148         static BOOST_UBLAS_INLINE
 149         result_type apply (argument1_type t1, argument2_type t2) {
 150             return t1 - t2;
 151         }
 152     };
 153     template<class T1, class T2>
 154     struct scalar_multiplies:
 155         public scalar_binary_functor<T1, T2> {
 156         typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;
 157         typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;
 158         typedef typename scalar_binary_functor<T1, T2>::result_type result_type;
 159
 160         static BOOST_UBLAS_INLINE
 161         result_type apply (argument1_type t1, argument2_type t2) {
 162             return t1 * t2;
 163         }
 164     };
 165     template<class T1, class T2>
 166     struct scalar_divides:
 167         public scalar_binary_functor<T1, T2> {
 168         typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type;
 169         typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type;
 170         typedef typename scalar_binary_functor<T1, T2>::result_type result_type;
 171
 172         static BOOST_UBLAS_INLINE
 173         result_type apply (argument1_type t1, argument2_type t2) {
 174             return t1 / t2;
 175         }
 176     };
 177
 178     template<class T1, class T2>
 179     struct scalar_binary_assign_functor {
 180         // ISSUE Remove reference to avoid reference to reference problems
 181         typedef typename type_traits<typename boost::remove_reference<T1>::type>::reference argument1_type;
 182         typedef typename type_traits<T2>::const_reference argument2_type;
 183     };
 184
 185     struct assign_tag {};
 186     struct computed_assign_tag {};
 187
 188     template<class T1, class T2>
 189     struct scalar_assign:
 190         public scalar_binary_assign_functor<T1, T2> {
 191         typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
 192         typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
 193 #if BOOST_WORKAROUND( __IBMCPP__, <=600 )
 194         static const bool computed ;
 195 #else
 196         static const bool computed = false ;
 197 #endif
 198
 199         static BOOST_UBLAS_INLINE
 200         void apply (argument1_type t1, argument2_type t2) {
 201             t1 = t2;
 202         }
 203
 204         template<class U1, class U2>
 205         struct rebind {
 206             typedef scalar_assign<U1, U2> other;
 207         };
 208     };
 209
 210 #if BOOST_WORKAROUND( __IBMCPP__, <=600 )
 211     template<class T1, class T2>
 212     const bool scalar_assign<T1,T2>::computed = false;
 213 #endif
 214
 215     template<class T1, class T2>
 216     struct scalar_plus_assign:
 217         public scalar_binary_assign_functor<T1, T2> {
 218         typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
 219         typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
 220 #if BOOST_WORKAROUND( __IBMCPP__, <=600 )
 221         static const bool computed ;
 222 #else
 223         static const bool computed = true ;
 224 #endif
 225
 226         static BOOST_UBLAS_INLINE
 227         void apply (argument1_type t1, argument2_type t2) {
 228             t1 += t2;
 229         }
 230
 231         template<class U1, class U2>
 232         struct rebind {
 233             typedef scalar_plus_assign<U1, U2> other;
 234         };
 235     };
 236
 237 #if BOOST_WORKAROUND( __IBMCPP__, <=600 )
 238     template<class T1, class T2>
 239     const bool scalar_plus_assign<T1,T2>::computed = true;
 240 #endif
 241
 242     template<class T1, class T2>
 243     struct scalar_minus_assign:
 244         public scalar_binary_assign_functor<T1, T2> {
 245         typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
 246         typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
 247 #if BOOST_WORKAROUND( __IBMCPP__, <=600 )
 248         static const bool computed ;
 249 #else
 250         static const bool computed = true ;
 251 #endif
 252
 253         static BOOST_UBLAS_INLINE
 254         void apply (argument1_type t1, argument2_type t2) {
 255             t1 -= t2;
 256         }
 257
 258         template<class U1, class U2>
 259         struct rebind {
 260             typedef scalar_minus_assign<U1, U2> other;
 261         };
 262     };
 263
 264 #if BOOST_WORKAROUND( __IBMCPP__, <=600 )
 265     template<class T1, class T2>
 266     const bool scalar_minus_assign<T1,T2>::computed = true;
 267 #endif
 268
 269     template<class T1, class T2>
 270     struct scalar_multiplies_assign:
 271         public scalar_binary_assign_functor<T1, T2> {
 272         typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
 273         typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
 274         static const bool computed = true;
 275
 276         static BOOST_UBLAS_INLINE
 277         void apply (argument1_type t1, argument2_type t2) {
 278             t1 *= t2;
 279         }
 280
 281         template<class U1, class U2>
 282         struct rebind {
 283             typedef scalar_multiplies_assign<U1, U2> other;
 284         };
 285     };
 286     template<class T1, class T2>
 287     struct scalar_divides_assign:
 288         public scalar_binary_assign_functor<T1, T2> {
 289         typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type;
 290         typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type;
 291         static const bool computed ;
 292
 293         static BOOST_UBLAS_INLINE
 294         void apply (argument1_type t1, argument2_type t2) {
 295             t1 /= t2;
 296         }
 297
 298         template<class U1, class U2>
 299         struct rebind {
 300             typedef scalar_divides_assign<U1, U2> other;
 301         };
 302     };
 303     template<class T1, class T2>
 304     const bool scalar_divides_assign<T1,T2>::computed = true;
 305
 306     template<class T1, class T2>
 307     struct scalar_binary_swap_functor {
 308         typedef typename type_traits<typename boost::remove_reference<T1>::type>::reference argument1_type;
 309         typedef typename type_traits<typename boost::remove_reference<T2>::type>::reference argument2_type;
 310     };
 311
 312     template<class T1, class T2>
 313     struct scalar_swap:
 314         public scalar_binary_swap_functor<T1, T2> {
 315         typedef typename scalar_binary_swap_functor<T1, T2>::argument1_type argument1_type;
 316         typedef typename scalar_binary_swap_functor<T1, T2>::argument2_type argument2_type;
 317
 318         static BOOST_UBLAS_INLINE
 319         void apply (argument1_type t1, argument2_type t2) {
 320             std::swap (t1, t2);
 321         }
 322
 323         template<class U1, class U2>
 324         struct rebind {
 325             typedef scalar_swap<U1, U2> other;
 326         };
 327     };
 328
 329     // Vector functors
 330
 331     // Unary returning scalar
 332     template<class V>
 333     struct vector_scalar_unary_functor {
 334         typedef typename V::value_type value_type;
 335         typedef typename V::value_type result_type;
 336     };
 337
 338     template<class V>
 339     struct vector_sum:
 340         public vector_scalar_unary_functor<V> {
 341         typedef typename vector_scalar_unary_functor<V>::value_type value_type;
 342         typedef typename vector_scalar_unary_functor<V>::result_type result_type;
 343
 344         template<class E>
 345         static BOOST_UBLAS_INLINE
 346         result_type apply (const vector_expression<E> &e) {
 347             result_type t = result_type (0);
 348             typedef typename E::size_type vector_size_type;
 349             vector_size_type size (e ().size ());
 350             for (vector_size_type i = 0; i < size; ++ i)
 351                 t += e () (i);
 352             return t;
 353         }
 354         // Dense case
 355         template<class D, class I>
 356         static BOOST_UBLAS_INLINE
 357         result_type apply (D size, I it) {
 358             result_type t = result_type (0);
 359             while (-- size >= 0)
 360                 t += *it, ++ it;
 361             return t;
 362         }
 363         // Sparse case
 364         template<class I>
 365         static BOOST_UBLAS_INLINE
 366         result_type apply (I it, const I &it_end) {
 367             result_type t = result_type (0);
 368             while (it != it_end)
 369                 t += *it, ++ it;
 370             return t;
 371         }
 372     };
 373
 374     // Unary returning real scalar
 375     template<class V>
 376     struct vector_scalar_real_unary_functor {
 377         typedef typename V::value_type value_type;
 378         typedef typename type_traits<value_type>::real_type real_type;
 379         typedef real_type result_type;
 380     };
 381
 382     template<class V>
 383     struct vector_norm_1:
 384         public vector_scalar_real_unary_functor<V> {
 385         typedef typename vector_scalar_real_unary_functor<V>::value_type value_type;
 386         typedef typename vector_scalar_real_unary_functor<V>::real_type real_type;
 387         typedef typename vector_scalar_real_unary_functor<V>::result_type result_type;
 388
 389         template<class E>
 390         static BOOST_UBLAS_INLINE
 391         result_type apply (const vector_expression<E> &e) {
 392             real_type t = real_type ();
 393             typedef typename E::size_type vector_size_type;
 394             vector_size_type size (e ().size ());
 395             for (vector_size_type i = 0; i < size; ++ i) {
 396                 real_type u (type_traits<value_type>::type_abs (e () (i)));
 397                 t += u;
 398             }
 399             return t;
 400         }
 401         // Dense case
 402         template<class D, class I>
 403         static BOOST_UBLAS_INLINE
 404         result_type apply (D size, I it) {
 405             real_type t = real_type ();
 406             while (-- size >= 0) {
 407                 real_type u (type_traits<value_type>::norm_1 (*it));
 408                 t += u;
 409                 ++ it;
 410             }
 411             return t;
 412         }
 413         // Sparse case
 414         template<class I>
 415         static BOOST_UBLAS_INLINE
 416         result_type apply (I it, const I &it_end) {
 417             real_type t = real_type ();
 418             while (it != it_end) {
 419                 real_type u (type_traits<value_type>::norm_1 (*it));
 420                 t += u;
 421                 ++ it;
 422             }
 423             return t;
 424         }
 425     };
 426     template<class V>
 427     struct vector_norm_2:
 428         public vector_scalar_real_unary_functor<V> {
 429         typedef typename vector_scalar_real_unary_functor<V>::value_type value_type;
 430         typedef typename vector_scalar_real_unary_functor<V>::real_type real_type;
 431         typedef typename vector_scalar_real_unary_functor<V>::result_type result_type;
 432
 433         template<class E>
 434         static BOOST_UBLAS_INLINE
 435         result_type apply (const vector_expression<E> &e) {
 436 #ifndef BOOST_UBLAS_SCALED_NORM
 437             real_type t = real_type ();
 438             typedef typename E::size_type vector_size_type;
 439             vector_size_type size (e ().size ());
 440             for (vector_size_type i = 0; i < size; ++ i) {
 441                 real_type u (type_traits<value_type>::norm_2 (e () (i)));
 442                 t +=  u * u;
 443             }
 444             return type_traits<real_type>::type_sqrt (t);
 445 #else
 446             real_type scale = real_type ();
 447             real_type sum_squares (1);
 448             size_type size (e ().size ());
 449             for (size_type i = 0; i < size; ++ i) {
 450                 real_type u (type_traits<value_type>::norm_2 (e () (i)));
 451                 if ( real_type () /* zero */ == u ) continue;
 452                 if (scale < u) {
 453                     real_type v (scale / u);
 454                     sum_squares = sum_squares * v * v + real_type (1);
 455                     scale = u;
 456                 } else {
 457                     real_type v (u / scale);
 458                     sum_squares += v * v;
 459                 }
 460             }
 461             return scale * type_traits<real_type>::type_sqrt (sum_squares);
 462 #endif
 463         }
 464         // Dense case
 465         template<class D, class I>
 466         static BOOST_UBLAS_INLINE
 467         result_type apply (D size, I it) {
 468 #ifndef BOOST_UBLAS_SCALED_NORM
 469             real_type t = real_type ();
 470             while (-- size >= 0) {
 471                 real_type u (type_traits<value_type>::norm_2 (*it));
 472                 t +=  u * u;
 473                 ++ it;
 474             }
 475             return type_traits<real_type>::type_sqrt (t);
 476 #else
 477             real_type scale = real_type ();
 478             real_type sum_squares (1);
 479             while (-- size >= 0) {
 480                 real_type u (type_traits<value_type>::norm_2 (*it));
 481                 if (scale < u) {
 482                     real_type v (scale / u);
 483                     sum_squares = sum_squares * v * v + real_type (1);
 484                     scale = u;
 485                 } else {
 486                     real_type v (u / scale);
 487                     sum_squares += v * v;
 488                 }
 489                 ++ it;
 490             }
 491             return scale * type_traits<real_type>::type_sqrt (sum_squares);
 492 #endif
 493         }
 494         // Sparse case
 495         template<class I>
 496         static BOOST_UBLAS_INLINE
 497         result_type apply (I it, const I &it_end) {
 498 #ifndef BOOST_UBLAS_SCALED_NORM
 499             real_type t = real_type ();
 500             while (it != it_end) {
 501                 real_type u (type_traits<value_type>::norm_2 (*it));
 502                 t +=  u * u;
 503                 ++ it;
 504             }
 505             return type_traits<real_type>::type_sqrt (t);
 506 #else
 507             real_type scale = real_type ();
 508             real_type sum_squares (1);
 509             while (it != it_end) {
 510                 real_type u (type_traits<value_type>::norm_2 (*it));
 511                 if (scale < u) {
 512                     real_type v (scale / u);
 513                     sum_squares = sum_squares * v * v + real_type (1);
 514                     scale = u;
 515                 } else {
 516                     real_type v (u / scale);
 517                     sum_squares += v * v;
 518                 }
 519                 ++ it;
 520             }
 521             return scale * type_traits<real_type>::type_sqrt (sum_squares);
 522 #endif
 523         }
 524     };
 525     template<class V>
 526     struct vector_norm_inf:
 527         public vector_scalar_real_unary_functor<V> {
 528         typedef typename vector_scalar_real_unary_functor<V>::value_type value_type;
 529         typedef typename vector_scalar_real_unary_functor<V>::real_type real_type;
 530         typedef typename vector_scalar_real_unary_functor<V>::result_type result_type;
 531
 532         template<class E>
 533         static BOOST_UBLAS_INLINE
 534         result_type apply (const vector_expression<E> &e) {
 535             real_type t = real_type ();
 536             typedef typename E::size_type vector_size_type;
 537             vector_size_type size (e ().size ());
 538             for (vector_size_type i = 0; i < size; ++ i) {
 539                 real_type u (type_traits<value_type>::norm_inf (e () (i)));
 540                 if (u > t)
 541                     t = u;
 542             }
 543             return t;
 544         }
 545         // Dense case
 546         template<class D, class I>
 547         static BOOST_UBLAS_INLINE
 548         result_type apply (D size, I it) {
 549             real_type t = real_type ();
 550             while (-- size >= 0) {
 551                 real_type u (type_traits<value_type>::norm_inf (*it));
 552                 if (u > t)
 553                     t = u;
 554                 ++ it;
 555             }
 556             return t;
 557         }
 558         // Sparse case
 559         template<class I>
 560         static BOOST_UBLAS_INLINE
 561         result_type apply (I it, const I &it_end) {
 562             real_type t = real_type ();
 563             while (it != it_end) {
 564                 real_type u (type_traits<value_type>::norm_inf (*it));
 565                 if (u > t)
 566                     t = u;
 567                 ++ it;
 568             }
 569             return t;
 570         }
 571     };
 572
 573     // Unary returning index
 574     template<class V>
 575     struct vector_scalar_index_unary_functor {
 576         typedef typename V::value_type value_type;
 577         typedef typename type_traits<value_type>::real_type real_type;
 578         typedef typename V::size_type result_type;
 579     };
 580
 581     template<class V>
 582     struct vector_index_norm_inf:
 583         public vector_scalar_index_unary_functor<V> {
 584         typedef typename vector_scalar_index_unary_functor<V>::value_type value_type;
 585         typedef typename vector_scalar_index_unary_functor<V>::real_type real_type;
 586         typedef typename vector_scalar_index_unary_functor<V>::result_type result_type;
 587
 588         template<class E>
 589         static BOOST_UBLAS_INLINE
 590         result_type apply (const vector_expression<E> &e) {
 591             // ISSUE For CBLAS compatibility return 0 index in empty case
 592             result_type i_norm_inf (0);
 593             real_type t = real_type ();
 594             typedef typename E::size_type vector_size_type;
 595             vector_size_type size (e ().size ());
 596             for (vector_size_type i = 0; i < size; ++ i) {
 597                 real_type u (type_traits<value_type>::norm_inf (e () (i)));
 598                 if (u > t) {
 599                     i_norm_inf = i;
 600                     t = u;
 601                 }
 602             }
 603             return i_norm_inf;
 604         }
 605         // Dense case
 606         template<class D, class I>
 607         static BOOST_UBLAS_INLINE
 608         result_type apply (D size, I it) {
 609             // ISSUE For CBLAS compatibility return 0 index in empty case
 610             result_type i_norm_inf (0);
 611             real_type t = real_type ();
 612             while (-- size >= 0) {
 613                 real_type u (type_traits<value_type>::norm_inf (*it));
 614                 if (u > t) {
 615                     i_norm_inf = it.index ();
 616                     t = u;
 617                 }
 618                 ++ it;
 619             }
 620             return i_norm_inf;
 621         }
 622         // Sparse case
 623         template<class I>
 624         static BOOST_UBLAS_INLINE
 625         result_type apply (I it, const I &it_end) {
 626             // ISSUE For CBLAS compatibility return 0 index in empty case
 627             result_type i_norm_inf (0);
 628             real_type t = real_type ();
 629             while (it != it_end) {
 630                 real_type u (type_traits<value_type>::norm_inf (*it));
 631                 if (u > t) {
 632                     i_norm_inf = it.index ();
 633                     t = u;
 634                 }
 635                 ++ it;
 636             }
 637             return i_norm_inf;
 638         }
 639     };
 640
 641     // Binary returning scalar
 642     template<class V1, class V2, class TV>
 643     struct vector_scalar_binary_functor {
 644         typedef TV value_type;
 645         typedef TV result_type;
 646     };
 647
 648     template<class V1, class V2, class TV>
 649     struct vector_inner_prod:
 650         public vector_scalar_binary_functor<V1, V2, TV> {
 651         typedef typename vector_scalar_binary_functor<V1, V2, TV>::value_type value_type;
 652         typedef typename vector_scalar_binary_functor<V1, V2, TV>::result_type result_type;
 653
 654         template<class C1, class C2>
 655         static BOOST_UBLAS_INLINE
 656         result_type apply (const vector_container<C1> &c1,
 657                            const vector_container<C2> &c2) {
 658 #ifdef BOOST_UBLAS_USE_SIMD
 659             using namespace raw;
 660             typedef typename C1::size_type vector_size_type;
 661             vector_size_type size (BOOST_UBLAS_SAME (c1 ().size (), c2 ().size ()));
 662             const typename V1::value_type *data1 = data_const (c1 ());
 663             const typename V1::value_type *data2 = data_const (c2 ());
 664             vector_size_type s1 = stride (c1 ());
 665             vector_size_type s2 = stride (c2 ());
 666             result_type t = result_type (0);
 667             if (s1 == 1 && s2 == 1) {
 668                 for (vector_size_type i = 0; i < size; ++ i)
 669                     t += data1 [i] * data2 [i];
 670             } else if (s2 == 1) {
 671                 for (vector_size_type i = 0, i1 = 0; i < size; ++ i, i1 += s1)
 672                     t += data1 [i1] * data2 [i];
 673             } else if (s1 == 1) {
 674                 for (vector_size_type i = 0, i2 = 0; i < size; ++ i, i2 += s2)
 675                     t += data1 [i] * data2 [i2];
 676             } else {
 677                 for (vector_size_type i = 0, i1 = 0, i2 = 0; i < size; ++ i, i1 += s1, i2 += s2)
 678                     t += data1 [i1] * data2 [i2];
 679             }
 680             return t;
 681 #elif defined(BOOST_UBLAS_HAVE_BINDINGS)
 682             return boost::numeric::bindings::atlas::dot (c1 (), c2 ());
 683 #else
 684             return apply (static_cast<const vector_expression<C1> > (c1), static_cast<const vector_expression<C2> > (c2));
 685 #endif
 686         }
 687         template<class E1, class E2>
 688         static BOOST_UBLAS_INLINE
 689         result_type apply (const vector_expression<E1> &e1,
 690                            const vector_expression<E2> &e2) {
 691             typedef typename E1::size_type vector_size_type;
 692             vector_size_type size (BOOST_UBLAS_SAME (e1 ().size (), e2 ().size ()));
 693             result_type t = result_type (0);
 694 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
 695             for (vector_size_type i = 0; i < size; ++ i)
 696                 t += e1 () (i) * e2 () (i);
 697 #else
 698             vector_size_type i (0);
 699             DD (size, 4, r, (t += e1 () (i) * e2 () (i), ++ i));
 700 #endif
 701             return t;
 702         }
 703         // Dense case
 704         template<class D, class I1, class I2>
 705         static BOOST_UBLAS_INLINE
 706         result_type apply (D size, I1 it1, I2 it2) {
 707             result_type t = result_type (0);
 708 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
 709             while (-- size >= 0)
 710                 t += *it1 * *it2, ++ it1, ++ it2;
 711 #else
 712             DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));
 713 #endif
 714             return t;
 715         }
 716         // Packed case
 717         template<class I1, class I2>
 718         static BOOST_UBLAS_INLINE
 719         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) {
 720             result_type t = result_type (0);
 721             typedef typename I1::difference_type vector_difference_type;
 722             vector_difference_type it1_size (it1_end - it1);
 723             vector_difference_type it2_size (it2_end - it2);
 724             vector_difference_type diff (0);
 725             if (it1_size > 0 && it2_size > 0)
 726                 diff = it2.index () - it1.index ();
 727             if (diff != 0) {
 728                 vector_difference_type size = (std::min) (diff, it1_size);
 729                 if (size > 0) {
 730                     it1 += size;
 731                     it1_size -= size;
 732                     diff -= size;
 733                 }
 734                 size = (std::min) (- diff, it2_size);
 735                 if (size > 0) {
 736                     it2 += size;
 737                     it2_size -= size;
 738                     diff += size;
 739                 }
 740             }
 741             vector_difference_type size ((std::min) (it1_size, it2_size));
 742             while (-- size >= 0)
 743                 t += *it1 * *it2, ++ it1, ++ it2;
 744             return t;
 745         }
 746         // Sparse case
 747         template<class I1, class I2>
 748         static BOOST_UBLAS_INLINE
 749         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) {
 750             result_type t = result_type (0);
 751             if (it1 != it1_end && it2 != it2_end) {
 752                 while (true) {
 753                     if (it1.index () == it2.index ()) {
 754                         t += *it1 * *it2, ++ it1, ++ it2;
 755                         if (it1 == it1_end || it2 == it2_end)
 756                             break;
 757                     } else if (it1.index () < it2.index ()) {
 758                         increment (it1, it1_end, it2.index () - it1.index ());
 759                         if (it1 == it1_end)
 760                             break;
 761                     } else if (it1.index () > it2.index ()) {
 762                         increment (it2, it2_end, it1.index () - it2.index ());
 763                         if (it2 == it2_end)
 764                             break;
 765                     }
 766                 }
 767             }
 768             return t;
 769         }
 770     };
 771
 772     // Matrix functors
 773
 774     // Binary returning vector
 775     template<class M1, class M2, class TV>
 776     struct matrix_vector_binary_functor {
 777         typedef typename M1::size_type size_type;
 778         typedef typename M1::difference_type difference_type;
 779         typedef TV value_type;
 780         typedef TV result_type;
 781     };
 782
 783     template<class M1, class M2, class TV>
 784     struct matrix_vector_prod1:
 785         public matrix_vector_binary_functor<M1, M2, TV> {
 786         typedef typename matrix_vector_binary_functor<M1, M2, TV>::size_type size_type;
 787         typedef typename matrix_vector_binary_functor<M1, M2, TV>::difference_type difference_type;
 788         typedef typename matrix_vector_binary_functor<M1, M2, TV>::value_type value_type;
 789         typedef typename matrix_vector_binary_functor<M1, M2, TV>::result_type result_type;
 790
 791         template<class C1, class C2>
 792         static BOOST_UBLAS_INLINE
 793         result_type apply (const matrix_container<C1> &c1,
 794                            const vector_container<C2> &c2,
 795                            size_type i) {
 796 #ifdef BOOST_UBLAS_USE_SIMD
 797             using namespace raw;
 798             size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().size ());
 799             const typename M1::value_type *data1 = data_const (c1 ()) + i * stride1 (c1 ());
 800             const typename M2::value_type *data2 = data_const (c2 ());
 801             size_type s1 = stride2 (c1 ());
 802             size_type s2 = stride (c2 ());
 803             result_type t = result_type (0);
 804             if (s1 == 1 && s2 == 1) {
 805                 for (size_type j = 0; j < size; ++ j)
 806                     t += data1 [j] * data2 [j];
 807             } else if (s2 == 1) {
 808                 for (size_type j = 0, j1 = 0; j < size; ++ j, j1 += s1)
 809                     t += data1 [j1] * data2 [j];
 810             } else if (s1 == 1) {
 811                 for (size_type j = 0, j2 = 0; j < size; ++ j, j2 += s2)
 812                     t += data1 [j] * data2 [j2];
 813             } else {
 814                 for (size_type j = 0, j1 = 0, j2 = 0; j < size; ++ j, j1 += s1, j2 += s2)
 815                     t += data1 [j1] * data2 [j2];
 816             }
 817             return t;
 818 #elif defined(BOOST_UBLAS_HAVE_BINDINGS)
 819             return boost::numeric::bindings::atlas::dot (c1 ().row (i), c2 ());
 820 #else
 821             return apply (static_cast<const matrix_expression<C1> > (c1), static_cast<const vector_expression<C2> > (c2, i));
 822 #endif
 823         }
 824         template<class E1, class E2>
 825         static BOOST_UBLAS_INLINE
 826         result_type apply (const matrix_expression<E1> &e1,
 827                            const vector_expression<E2> &e2,
 828                            size_type i) {
 829             size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size ());
 830             result_type t = result_type (0);
 831 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
 832             for (size_type j = 0; j < size; ++ j)
 833                 t += e1 () (i, j) * e2 () (j);
 834 #else
 835             size_type j (0);
 836             DD (size, 4, r, (t += e1 () (i, j) * e2 () (j), ++ j));
 837 #endif
 838             return t;
 839         }
 840         // Dense case
 841         template<class I1, class I2>
 842         static BOOST_UBLAS_INLINE
 843         result_type apply (difference_type size, I1 it1, I2 it2) {
 844             result_type t = result_type (0);
 845 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
 846             while (-- size >= 0)
 847                 t += *it1 * *it2, ++ it1, ++ it2;
 848 #else
 849             DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));
 850 #endif
 851             return t;
 852         }
 853         // Packed case
 854         template<class I1, class I2>
 855         static BOOST_UBLAS_INLINE
 856         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) {
 857             result_type t = result_type (0);
 858             difference_type it1_size (it1_end - it1);
 859             difference_type it2_size (it2_end - it2);
 860             difference_type diff (0);
 861             if (it1_size > 0 && it2_size > 0)
 862                 diff = it2.index () - it1.index2 ();
 863             if (diff != 0) {
 864                 difference_type size = (std::min) (diff, it1_size);
 865                 if (size > 0) {
 866                     it1 += size;
 867                     it1_size -= size;
 868                     diff -= size;
 869                 }
 870                 size = (std::min) (- diff, it2_size);
 871                 if (size > 0) {
 872                     it2 += size;
 873                     it2_size -= size;
 874                     diff += size;
 875                 }
 876             }
 877             difference_type size ((std::min) (it1_size, it2_size));
 878             while (-- size >= 0)
 879                 t += *it1 * *it2, ++ it1, ++ it2;
 880             return t;
 881         }
 882         // Sparse case
 883         template<class I1, class I2>
 884         static BOOST_UBLAS_INLINE
 885         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
 886                            sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) {
 887             result_type t = result_type (0);
 888             if (it1 != it1_end && it2 != it2_end) {
 889                 size_type it1_index = it1.index2 (), it2_index = it2.index ();
 890                 while (true) {
 891                     difference_type compare = it1_index - it2_index;
 892                     if (compare == 0) {
 893                         t += *it1 * *it2, ++ it1, ++ it2;
 894                         if (it1 != it1_end && it2 != it2_end) {
 895                             it1_index = it1.index2 ();
 896                             it2_index = it2.index ();
 897                         } else
 898                             break;
 899                     } else if (compare < 0) {
 900                         increment (it1, it1_end, - compare);
 901                         if (it1 != it1_end)
 902                             it1_index = it1.index2 ();
 903                         else
 904                             break;
 905                     } else if (compare > 0) {
 906                         increment (it2, it2_end, compare);
 907                         if (it2 != it2_end)
 908                             it2_index = it2.index ();
 909                         else
 910                             break;
 911                     }
 912                 }
 913             }
 914             return t;
 915         }
 916         // Sparse packed case
 917         template<class I1, class I2>
 918         static BOOST_UBLAS_INLINE
 919         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */,
 920                            sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) {
 921             result_type t = result_type (0);
 922             while (it1 != it1_end) {
 923                 t += *it1 * it2 () (it1.index2 ());
 924                 ++ it1;
 925             }
 926             return t;
 927         }
 928         // Packed sparse case
 929         template<class I1, class I2>
 930         static BOOST_UBLAS_INLINE
 931         result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end,
 932                            packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) {
 933             result_type t = result_type (0);
 934             while (it2 != it2_end) {
 935                 t += it1 () (it1.index1 (), it2.index ()) * *it2;
 936                 ++ it2;
 937             }
 938             return t;
 939         }
 940         // Another dispatcher
 941         template<class I1, class I2>
 942         static BOOST_UBLAS_INLINE
 943         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
 944                            sparse_bidirectional_iterator_tag) {
 945             typedef typename I1::iterator_category iterator1_category;
 946             typedef typename I2::iterator_category iterator2_category;
 947             return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ());
 948         }
 949     };
 950
 951     template<class M1, class M2, class TV>
 952     struct matrix_vector_prod2:
 953         public matrix_vector_binary_functor<M1, M2, TV> {
 954         typedef typename matrix_vector_binary_functor<M1, M2, TV>::size_type size_type;
 955         typedef typename matrix_vector_binary_functor<M1, M2, TV>::difference_type difference_type;
 956         typedef typename matrix_vector_binary_functor<M1, M2, TV>::value_type value_type;
 957         typedef typename matrix_vector_binary_functor<M1, M2, TV>::result_type result_type;
 958
 959         template<class C1, class C2>
 960         static BOOST_UBLAS_INLINE
 961         result_type apply (const vector_container<C1> &c1,
 962                            const matrix_container<C2> &c2,
 963                            size_type i) {
 964 #ifdef BOOST_UBLAS_USE_SIMD
 965             using namespace raw;
 966             size_type size = BOOST_UBLAS_SAME (c1 ().size (), c2 ().size1 ());
 967             const typename M1::value_type *data1 = data_const (c1 ());
 968             const typename M2::value_type *data2 = data_const (c2 ()) + i * stride2 (c2 ());
 969             size_type s1 = stride (c1 ());
 970             size_type s2 = stride1 (c2 ());
 971             result_type t = result_type (0);
 972             if (s1 == 1 && s2 == 1) {
 973                 for (size_type j = 0; j < size; ++ j)
 974                     t += data1 [j] * data2 [j];
 975             } else if (s2 == 1) {
 976                 for (size_type j = 0, j1 = 0; j < size; ++ j, j1 += s1)
 977                     t += data1 [j1] * data2 [j];
 978             } else if (s1 == 1) {
 979                 for (size_type j = 0, j2 = 0; j < size; ++ j, j2 += s2)
 980                     t += data1 [j] * data2 [j2];
 981             } else {
 982                 for (size_type j = 0, j1 = 0, j2 = 0; j < size; ++ j, j1 += s1, j2 += s2)
 983                     t += data1 [j1] * data2 [j2];
 984             }
 985             return t;
 986 #elif defined(BOOST_UBLAS_HAVE_BINDINGS)
 987             return boost::numeric::bindings::atlas::dot (c1 (), c2 ().column (i));
 988 #else
 989             return apply (static_cast<const vector_expression<C1> > (c1), static_cast<const matrix_expression<C2> > (c2, i));
 990 #endif
 991         }
 992         template<class E1, class E2>
 993         static BOOST_UBLAS_INLINE
 994         result_type apply (const vector_expression<E1> &e1,
 995                            const matrix_expression<E2> &e2,
 996                            size_type i) {
 997             size_type size = BOOST_UBLAS_SAME (e1 ().size (), e2 ().size1 ());
 998             result_type t = result_type (0);
 999 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
1000             for (size_type j = 0; j < size; ++ j)
1001                 t += e1 () (j) * e2 () (j, i);
1002 #else
1003             size_type j (0);
1004             DD (size, 4, r, (t += e1 () (j) * e2 () (j, i), ++ j));
1005 #endif
1006             return t;
1007         }
1008         // Dense case
1009         template<class I1, class I2>
1010         static BOOST_UBLAS_INLINE
1011         result_type apply (difference_type size, I1 it1, I2 it2) {
1012             result_type t = result_type (0);
1013 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
1014             while (-- size >= 0)
1015                 t += *it1 * *it2, ++ it1, ++ it2;
1016 #else
1017             DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));
1018 #endif
1019             return t;
1020         }
1021         // Packed case
1022         template<class I1, class I2>
1023         static BOOST_UBLAS_INLINE
1024         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) {
1025             result_type t = result_type (0);
1026             difference_type it1_size (it1_end - it1);
1027             difference_type it2_size (it2_end - it2);
1028             difference_type diff (0);
1029             if (it1_size > 0 && it2_size > 0)
1030                 diff = it2.index1 () - it1.index ();
1031             if (diff != 0) {
1032                 difference_type size = (std::min) (diff, it1_size);
1033                 if (size > 0) {
1034                     it1 += size;
1035                     it1_size -= size;
1036                     diff -= size;
1037                 }
1038                 size = (std::min) (- diff, it2_size);
1039                 if (size > 0) {
1040                     it2 += size;
1041                     it2_size -= size;
1042                     diff += size;
1043                 }
1044             }
1045             difference_type size ((std::min) (it1_size, it2_size));
1046             while (-- size >= 0)
1047                 t += *it1 * *it2, ++ it1, ++ it2;
1048             return t;
1049         }
1050         // Sparse case
1051         template<class I1, class I2>
1052         static BOOST_UBLAS_INLINE
1053         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
1054                            sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) {
1055             result_type t = result_type (0);
1056             if (it1 != it1_end && it2 != it2_end) {
1057                 size_type it1_index = it1.index (), it2_index = it2.index1 ();
1058                 while (true) {
1059                     difference_type compare = it1_index - it2_index;
1060                     if (compare == 0) {
1061                         t += *it1 * *it2, ++ it1, ++ it2;
1062                         if (it1 != it1_end && it2 != it2_end) {
1063                             it1_index = it1.index ();
1064                             it2_index = it2.index1 ();
1065                         } else
1066                             break;
1067                     } else if (compare < 0) {
1068                         increment (it1, it1_end, - compare);
1069                         if (it1 != it1_end)
1070                             it1_index = it1.index ();
1071                         else
1072                             break;
1073                     } else if (compare > 0) {
1074                         increment (it2, it2_end, compare);
1075                         if (it2 != it2_end)
1076                             it2_index = it2.index1 ();
1077                         else
1078                             break;
1079                     }
1080                 }
1081             }
1082             return t;
1083         }
1084         // Packed sparse case
1085         template<class I1, class I2>
1086         static BOOST_UBLAS_INLINE
1087         result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end,
1088                            packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) {
1089             result_type t = result_type (0);
1090             while (it2 != it2_end) {
1091                 t += it1 () (it2.index1 ()) * *it2;
1092                 ++ it2;
1093             }
1094             return t;
1095         }
1096         // Sparse packed case
1097         template<class I1, class I2>
1098         static BOOST_UBLAS_INLINE
1099         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */,
1100                            sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) {
1101             result_type t = result_type (0);
1102             while (it1 != it1_end) {
1103                 t += *it1 * it2 () (it1.index (), it2.index2 ());
1104                 ++ it1;
1105             }
1106             return t;
1107         }
1108         // Another dispatcher
1109         template<class I1, class I2>
1110         static BOOST_UBLAS_INLINE
1111         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end,
1112                            sparse_bidirectional_iterator_tag) {
1113             typedef typename I1::iterator_category iterator1_category;
1114             typedef typename I2::iterator_category iterator2_category;
1115             return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ());
1116         }
1117     };
1118
1119     // Binary returning matrix
1120     template<class M1, class M2, class TV>
1121     struct matrix_matrix_binary_functor {
1122         typedef typename M1::size_type size_type;
1123         typedef typename M1::difference_type difference_type;
1124         typedef TV value_type;
1125         typedef TV result_type;
1126     };
1127
1128     template<class M1, class M2, class TV>
1129     struct matrix_matrix_prod:
1130         public matrix_matrix_binary_functor<M1, M2, TV> {
1131         typedef typename matrix_matrix_binary_functor<M1, M2, TV>::size_type size_type;
1132         typedef typename matrix_matrix_binary_functor<M1, M2, TV>::difference_type difference_type;
1133         typedef typename matrix_matrix_binary_functor<M1, M2, TV>::value_type value_type;
1134         typedef typename matrix_matrix_binary_functor<M1, M2, TV>::result_type result_type;
1135
1136         template<class C1, class C2>
1137         static BOOST_UBLAS_INLINE
1138         result_type apply (const matrix_container<C1> &c1,
1139                            const matrix_container<C2> &c2,
1140                            size_type i, size_type j) {
1141 #ifdef BOOST_UBLAS_USE_SIMD
1142             using namespace raw;
1143             size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().sizc1 ());
1144             const typename M1::value_type *data1 = data_const (c1 ()) + i * stride1 (c1 ());
1145             const typename M2::value_type *data2 = data_const (c2 ()) + j * stride2 (c2 ());
1146             size_type s1 = stride2 (c1 ());
1147             size_type s2 = stride1 (c2 ());
1148             result_type t = result_type (0);
1149             if (s1 == 1 && s2 == 1) {
1150                 for (size_type k = 0; k < size; ++ k)
1151                     t += data1 [k] * data2 [k];
1152             } else if (s2 == 1) {
1153                 for (size_type k = 0, k1 = 0; k < size; ++ k, k1 += s1)
1154                     t += data1 [k1] * data2 [k];
1155             } else if (s1 == 1) {
1156                 for (size_type k = 0, k2 = 0; k < size; ++ k, k2 += s2)
1157                     t += data1 [k] * data2 [k2];
1158             } else {
1159                 for (size_type k = 0, k1 = 0, k2 = 0; k < size; ++ k, k1 += s1, k2 += s2)
1160                     t += data1 [k1] * data2 [k2];
1161             }
1162             return t;
1163 #elif defined(BOOST_UBLAS_HAVE_BINDINGS)
1164             return boost::numeric::bindings::atlas::dot (c1 ().row (i), c2 ().column (j));
1165 #else
1166             boost::ignore_unused(j);
1167             return apply (static_cast<const matrix_expression<C1> > (c1), static_cast<const matrix_expression<C2> > (c2, i));
1168 #endif
1169         }
1170         template<class E1, class E2>
1171         static BOOST_UBLAS_INLINE
1172         result_type apply (const matrix_expression<E1> &e1,
1173                            const matrix_expression<E2> &e2,
1174                            size_type i, size_type j) {
1175             size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ());
1176             result_type t = result_type (0);
1177 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
1178             for (size_type k = 0; k < size; ++ k)
1179                 t += e1 () (i, k) * e2 () (k, j);
1180 #else
1181             size_type k (0);
1182             DD (size, 4, r, (t += e1 () (i, k) * e2 () (k, j), ++ k));
1183 #endif
1184             return t;
1185         }
1186         // Dense case
1187         template<class I1, class I2>
1188         static BOOST_UBLAS_INLINE
1189         result_type apply (difference_type size, I1 it1, I2 it2) {
1190             result_type t = result_type (0);
1191 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE
1192             while (-- size >= 0)
1193                 t += *it1 * *it2, ++ it1, ++ it2;
1194 #else
1195             DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2));
1196 #endif
1197             return t;
1198         }
1199         // Packed case
1200         template<class I1, class I2>
1201         static BOOST_UBLAS_INLINE
1202         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, packed_random_access_iterator_tag) {
1203             result_type t = result_type (0);
1204             difference_type it1_size (it1_end - it1);
1205             difference_type it2_size (it2_end - it2);
1206             difference_type diff (0);
1207             if (it1_size > 0 && it2_size > 0)
1208                 diff = it2.index1 () - it1.index2 ();
1209             if (diff != 0) {
1210                 difference_type size = (std::min) (diff, it1_size);
1211                 if (size > 0) {
1212                     it1 += size;
1213                     it1_size -= size;
1214                     diff -= size;
1215                 }
1216                 size = (std::min) (- diff, it2_size);
1217                 if (size > 0) {
1218                     it2 += size;
1219                     it2_size -= size;
1220                     diff += size;
1221                 }
1222             }
1223             difference_type size ((std::min) (it1_size, it2_size));
1224             while (-- size >= 0)
1225                 t += *it1 * *it2, ++ it1, ++ it2;
1226             return t;
1227         }
1228         // Sparse case
1229         template<class I1, class I2>
1230         static BOOST_UBLAS_INLINE
1231         result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) {
1232             result_type t = result_type (0);
1233             if (it1 != it1_end && it2 != it2_end) {
1234                 size_type it1_index = it1.index2 (), it2_index = it2.index1 ();
1235                 while (true) {
1236                     difference_type compare = difference_type (it1_index - it2_index);
1237                     if (compare == 0) {
1238                         t += *it1 * *it2, ++ it1, ++ it2;
1239                         if (it1 != it1_end && it2 != it2_end) {
1240                             it1_index = it1.index2 ();
1241                             it2_index = it2.index1 ();
1242                         } else
1243                             break;
1244                     } else if (compare < 0) {
1245                         increment (it1, it1_end, - compare);
1246                         if (it1 != it1_end)
1247                             it1_index = it1.index2 ();
1248                         else
1249                             break;
1250                     } else if (compare > 0) {
1251                         increment (it2, it2_end, compare);
1252                         if (it2 != it2_end)
1253                             it2_index = it2.index1 ();
1254                         else
1255                             break;
1256                     }
1257                 }
1258             }
1259             return t;
1260         }
1261     };
1262
1263     // Unary returning scalar norm
1264     template<class M>
1265     struct matrix_scalar_real_unary_functor {
1266         typedef typename M::value_type value_type;
1267         typedef typename type_traits<value_type>::real_type real_type;
1268         typedef real_type result_type;
1269     };
1270
1271     template<class M>
1272     struct matrix_norm_1:
1273         public matrix_scalar_real_unary_functor<M> {
1274         typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type;
1275         typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type;
1276         typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type;
1277
1278         template<class E>
1279         static BOOST_UBLAS_INLINE
1280         result_type apply (const matrix_expression<E> &e) {
1281             real_type t = real_type ();
1282             typedef typename E::size_type matrix_size_type;
1283             matrix_size_type size2 (e ().size2 ());
1284             for (matrix_size_type j = 0; j < size2; ++ j) {
1285                 real_type u = real_type ();
1286                 matrix_size_type size1 (e ().size1 ());
1287                 for (matrix_size_type i = 0; i < size1; ++ i) {
1288                     real_type v (type_traits<value_type>::norm_1 (e () (i, j)));
1289                     u += v;
1290                 }
1291                 if (u > t)
1292                     t = u;
1293             }
1294             return t;
1295         }
1296     };
1297
1298     template<class M>
1299     struct matrix_norm_frobenius:
1300         public matrix_scalar_real_unary_functor<M> {
1301         typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type;
1302         typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type;
1303         typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type;
1304
1305         template<class E>
1306         static BOOST_UBLAS_INLINE
1307         result_type apply (const matrix_expression<E> &e) {
1308             real_type t = real_type ();
1309             typedef typename E::size_type matrix_size_type;
1310             matrix_size_type size1 (e ().size1 ());
1311             for (matrix_size_type i = 0; i < size1; ++ i) {
1312                 matrix_size_type size2 (e ().size2 ());
1313                 for (matrix_size_type j = 0; j < size2; ++ j) {
1314                     real_type u (type_traits<value_type>::norm_2 (e () (i, j)));
1315                     t +=  u * u;
1316                 }
1317             }
1318             return type_traits<real_type>::type_sqrt (t);
1319         }
1320     };
1321
1322     template<class M>
1323     struct matrix_norm_inf:
1324         public matrix_scalar_real_unary_functor<M> {
1325         typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type;
1326         typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type;
1327         typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type;
1328
1329         template<class E>
1330         static BOOST_UBLAS_INLINE
1331         result_type apply (const matrix_expression<E> &e) {
1332             real_type t = real_type ();
1333             typedef typename E::size_type matrix_size_type;
1334             matrix_size_type size1 (e ().size1 ());
1335             for (matrix_size_type i = 0; i < size1; ++ i) {
1336                 real_type u = real_type ();
1337                 matrix_size_type size2 (e ().size2 ());
1338                 for (matrix_size_type j = 0; j < size2; ++ j) {
1339                     real_type v (type_traits<value_type>::norm_inf (e () (i, j)));
1340                     u += v;
1341                 }
1342                 if (u > t)
1343                     t = u;
1344             }
1345             return t;
1346         }
1347     };
1348
1349     // forward declaration
1350     template <class Z, class D> struct basic_column_major;
1351
1352     // This functor defines storage layout and it's properties
1353     // matrix (i,j) -> storage [i * size_i + j]
1354     template <class Z, class D>
1355     struct basic_row_major {
1356         typedef Z size_type;
1357         typedef D difference_type;
1358         typedef row_major_tag orientation_category;
1359         typedef basic_column_major<Z,D> transposed_layout;
1360
1361         static
1362         BOOST_UBLAS_INLINE
1363         size_type storage_size (size_type size_i, size_type size_j) {
1364             // Guard against size_type overflow
1365             BOOST_UBLAS_CHECK (size_j == 0 || size_i <= (std::numeric_limits<size_type>::max) () / size_j, bad_size ());
1366             return size_i * size_j;
1367         }
1368
1369         // Indexing conversion to storage element
1370         static
1371         BOOST_UBLAS_INLINE
1372         size_type element (size_type i, size_type size_i, size_type j, size_type size_j) {
1373             BOOST_UBLAS_CHECK (i < size_i, bad_index ());
1374             BOOST_UBLAS_CHECK (j < size_j, bad_index ());
1375             detail::ignore_unused_variable_warning(size_i);
1376             // Guard against size_type overflow
1377             BOOST_UBLAS_CHECK (i <= ((std::numeric_limits<size_type>::max) () - j) / size_j, bad_index ());
1378             return i * size_j + j;
1379         }
1380         static
1381         BOOST_UBLAS_INLINE
1382         size_type address (size_type i, size_type size_i, size_type j, size_type size_j) {
1383             BOOST_UBLAS_CHECK (i <= size_i, bad_index ());
1384             BOOST_UBLAS_CHECK (j <= size_j, bad_index ());
1385             // Guard against size_type overflow - address may be size_j past end of storage
1386             BOOST_UBLAS_CHECK (size_j == 0 || i <= ((std::numeric_limits<size_type>::max) () - j) / size_j, bad_index ());
1387             detail::ignore_unused_variable_warning(size_i);
1388             return i * size_j + j;
1389         }
1390
1391         // Storage element to index conversion
1392         static
1393         BOOST_UBLAS_INLINE
1394         difference_type distance_i (difference_type k, size_type /* size_i */, size_type size_j) {
1395             return size_j != 0 ? k / size_j : 0;
1396         }
1397         static
1398         BOOST_UBLAS_INLINE
1399         difference_type distance_j (difference_type k, size_type /* size_i */, size_type /* size_j */) {
1400             return k;
1401         }
1402         static
1403         BOOST_UBLAS_INLINE
1404         size_type index_i (difference_type k, size_type /* size_i */, size_type size_j) {
1405             return size_j != 0 ? k / size_j : 0;
1406         }
1407         static
1408         BOOST_UBLAS_INLINE
1409         size_type index_j (difference_type k, size_type /* size_i */, size_type size_j) {
1410             return size_j != 0 ? k % size_j : 0;
1411         }
1412         static
1413         BOOST_UBLAS_INLINE
1414         bool fast_i () {
1415             return false;
1416         }
1417         static
1418         BOOST_UBLAS_INLINE
1419         bool fast_j () {
1420             return true;
1421         }
1422
1423         // Iterating storage elements
1424         template<class I>
1425         static
1426         BOOST_UBLAS_INLINE
1427         void increment_i (I &it, size_type /* size_i */, size_type size_j) {
1428             it += size_j;
1429         }
1430         template<class I>
1431         static
1432         BOOST_UBLAS_INLINE
1433         void increment_i (I &it, difference_type n, size_type /* size_i */, size_type size_j) {
1434             it += n * size_j;
1435         }
1436         template<class I>
1437         static
1438         BOOST_UBLAS_INLINE
1439         void decrement_i (I &it, size_type /* size_i */, size_type size_j) {
1440             it -= size_j;
1441         }
1442         template<class I>
1443         static
1444         BOOST_UBLAS_INLINE
1445         void decrement_i (I &it, difference_type n, size_type /* size_i */, size_type size_j) {
1446             it -= n * size_j;
1447         }
1448         template<class I>
1449         static
1450         BOOST_UBLAS_INLINE
1451         void increment_j (I &it, size_type /* size_i */, size_type /* size_j */) {
1452             ++ it;
1453         }
1454         template<class I>
1455         static
1456         BOOST_UBLAS_INLINE
1457         void increment_j (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) {
1458             it += n;
1459         }
1460         template<class I>
1461         static
1462         BOOST_UBLAS_INLINE
1463         void decrement_j (I &it, size_type /* size_i */, size_type /* size_j */) {
1464             -- it;
1465         }
1466         template<class I>
1467         static
1468         BOOST_UBLAS_INLINE
1469         void decrement_j (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) {
1470             it -= n;
1471         }
1472
1473         // Triangular access
1474         static
1475         BOOST_UBLAS_INLINE
1476         size_type triangular_size (size_type size_i, size_type size_j) {
1477             size_type size = (std::max) (size_i, size_j);
1478             // Guard against size_type overflow - simplified
1479             BOOST_UBLAS_CHECK (size == 0 || size / 2 < (std::numeric_limits<size_type>::max) () / size /* +1/2 */, bad_size ());
1480             return ((size + 1) * size) / 2;
1481         }
1482         static
1483         BOOST_UBLAS_INLINE
1484         size_type lower_element (size_type i, size_type size_i, size_type j, size_type size_j) {
1485             BOOST_UBLAS_CHECK (i < size_i, bad_index ());
1486             BOOST_UBLAS_CHECK (j < size_j, bad_index ());
1487             BOOST_UBLAS_CHECK (i >= j, bad_index ());
1488             detail::ignore_unused_variable_warning(size_i);
1489             detail::ignore_unused_variable_warning(size_j);
1490             // FIXME size_type overflow
1491             // sigma_i (i + 1) = (i + 1) * i / 2
1492             // i = 0 1 2 3, sigma = 0 1 3 6
1493             return ((i + 1) * i) / 2 + j;
1494         }
1495         static
1496         BOOST_UBLAS_INLINE
1497         size_type upper_element (size_type i, size_type size_i, size_type j, size_type size_j) {
1498             BOOST_UBLAS_CHECK (i < size_i, bad_index ());
1499             BOOST_UBLAS_CHECK (j < size_j, bad_index ());
1500             BOOST_UBLAS_CHECK (i <= j, bad_index ());
1501             // FIXME size_type overflow
1502             // sigma_i (size - i) = size * i - i * (i - 1) / 2
1503             // i = 0 1 2 3, sigma = 0 4 7 9
1504             return (i * (2 * (std::max) (size_i, size_j) - i + 1)) / 2 + j - i;
1505         }
1506
1507         // Major and minor indices
1508         static
1509         BOOST_UBLAS_INLINE
1510         size_type index_M (size_type index1, size_type /* index2 */) {
1511             return index1;
1512         }
1513         static
1514         BOOST_UBLAS_INLINE
1515         size_type index_m (size_type /* index1 */, size_type index2) {
1516             return index2;
1517         }
1518         static
1519         BOOST_UBLAS_INLINE
1520         size_type size_M (size_type size_i, size_type /* size_j */) {
1521             return size_i;
1522         }
1523         static
1524         BOOST_UBLAS_INLINE
1525         size_type size_m (size_type /* size_i */, size_type size_j) {
1526             return size_j;
1527         }
1528     };
1529
1530     // This functor defines storage layout and it's properties
1531     // matrix (i,j) -> storage [i + j * size_i]
1532     template <class Z, class D>
1533     struct basic_column_major {
1534         typedef Z size_type;
1535         typedef D difference_type;
1536         typedef column_major_tag orientation_category;
1537         typedef basic_row_major<Z,D> transposed_layout;
1538
1539         static
1540         BOOST_UBLAS_INLINE
1541         size_type storage_size (size_type size_i, size_type size_j) {
1542             // Guard against size_type overflow
1543             BOOST_UBLAS_CHECK (size_i == 0 || size_j <= (std::numeric_limits<size_type>::max) () / size_i, bad_size ());
1544             return size_i * size_j;
1545         }
1546
1547         // Indexing conversion to storage element
1548         static
1549         BOOST_UBLAS_INLINE
1550         size_type element (size_type i, size_type size_i, size_type j, size_type size_j) {
1551             BOOST_UBLAS_CHECK (i < size_i, bad_index ());
1552             BOOST_UBLAS_CHECK (j < size_j, bad_index ());
1553             detail::ignore_unused_variable_warning(size_j);
1554             // Guard against size_type overflow
1555             BOOST_UBLAS_CHECK (j <= ((std::numeric_limits<size_type>::max) () - i) / size_i, bad_index ());
1556             return i + j * size_i;
1557         }
1558         static
1559         BOOST_UBLAS_INLINE
1560         size_type address (size_type i, size_type size_i, size_type j, size_type size_j) {
1561             BOOST_UBLAS_CHECK (i <= size_i, bad_index ());
1562             BOOST_UBLAS_CHECK (j <= size_j, bad_index ());
1563             detail::ignore_unused_variable_warning(size_j);
1564             // Guard against size_type overflow - address may be size_i past end of storage
1565             BOOST_UBLAS_CHECK (size_i == 0 || j <= ((std::numeric_limits<size_type>::max) () - i) / size_i, bad_index ());
1566             return i + j * size_i;
1567         }
1568
1569         // Storage element to index conversion
1570         static
1571         BOOST_UBLAS_INLINE
1572         difference_type distance_i (difference_type k, size_type /* size_i */, size_type /* size_j */) {
1573             return k;
1574         }
1575         static
1576         BOOST_UBLAS_INLINE
1577         difference_type distance_j (difference_type k, size_type size_i, size_type /* size_j */) {
1578             return size_i != 0 ? k / size_i : 0;
1579         }
1580         static
1581         BOOST_UBLAS_INLINE
1582         size_type index_i (difference_type k, size_type size_i, size_type /* size_j */) {
1583             return size_i != 0 ? k % size_i : 0;
1584         }
1585         static
1586         BOOST_UBLAS_INLINE
1587         size_type index_j (difference_type k, size_type size_i, size_type /* size_j */) {
1588             return size_i != 0 ? k / size_i : 0;
1589         }
1590         static
1591         BOOST_UBLAS_INLINE
1592         bool fast_i () {
1593             return true;
1594         }
1595         static
1596         BOOST_UBLAS_INLINE
1597         bool fast_j () {
1598             return false;
1599         }
1600
1601         // Iterating
1602         template<class I>
1603         static
1604         BOOST_UBLAS_INLINE
1605         void increment_i (I &it, size_type /* size_i */, size_type /* size_j */) {
1606             ++ it;
1607         }
1608         template<class I>
1609         static
1610         BOOST_UBLAS_INLINE
1611         void increment_i (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) {
1612             it += n;
1613         }
1614         template<class I>
1615         static
1616         BOOST_UBLAS_INLINE
1617         void decrement_i (I &it, size_type /* size_i */, size_type /* size_j */) {
1618             -- it;
1619         }
1620         template<class I>
1621         static
1622         BOOST_UBLAS_INLINE
1623         void decrement_i (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) {
1624             it -= n;
1625         }
1626         template<class I>
1627         static
1628         BOOST_UBLAS_INLINE
1629         void increment_j (I &it, size_type size_i, size_type /* size_j */) {
1630             it += size_i;
1631         }
1632         template<class I>
1633         static
1634         BOOST_UBLAS_INLINE
1635         void increment_j (I &it, difference_type n, size_type size_i, size_type /* size_j */) {
1636             it += n * size_i;
1637         }
1638         template<class I>
1639         static
1640         BOOST_UBLAS_INLINE
1641         void decrement_j (I &it, size_type size_i, size_type /* size_j */) {
1642             it -= size_i;
1643         }
1644         template<class I>
1645         static
1646         BOOST_UBLAS_INLINE
1647         void decrement_j (I &it, difference_type n, size_type size_i, size_type /* size_j */) {
1648             it -= n* size_i;
1649         }
1650
1651         // Triangular access
1652         static
1653         BOOST_UBLAS_INLINE
1654         size_type triangular_size (size_type size_i, size_type size_j) {
1655             size_type size = (std::max) (size_i, size_j);
1656             // Guard against size_type overflow - simplified
1657             BOOST_UBLAS_CHECK (size == 0 || size / 2 < (std::numeric_limits<size_type>::max) () / size /* +1/2 */, bad_size ());
1658             return ((size + 1) * size) / 2;
1659         }
1660         static
1661         BOOST_UBLAS_INLINE
1662         size_type lower_element (size_type i, size_type size_i, size_type j, size_type size_j) {
1663             BOOST_UBLAS_CHECK (i < size_i, bad_index ());
1664             BOOST_UBLAS_CHECK (j < size_j, bad_index ());
1665             BOOST_UBLAS_CHECK (i >= j, bad_index ());
1666             // FIXME size_type overflow
1667             // sigma_j (size - j) = size * j - j * (j - 1) / 2
1668             // j = 0 1 2 3, sigma = 0 4 7 9
1669             return i - j + (j * (2 * (std::max) (size_i, size_j) - j + 1)) / 2;
1670         }
1671         static
1672         BOOST_UBLAS_INLINE
1673         size_type upper_element (size_type i, size_type size_i, size_type j, size_type size_j) {
1674             BOOST_UBLAS_CHECK (i < size_i, bad_index ());
1675             BOOST_UBLAS_CHECK (j < size_j, bad_index ());
1676             BOOST_UBLAS_CHECK (i <= j, bad_index ());
1677             // FIXME size_type overflow
1678             // sigma_j (j + 1) = (j + 1) * j / 2
1679             // j = 0 1 2 3, sigma = 0 1 3 6
1680             return i + ((j + 1) * j) / 2;
1681         }
1682
1683         // Major and minor indices
1684         static
1685         BOOST_UBLAS_INLINE
1686         size_type index_M (size_type /* index1 */, size_type index2) {
1687             return index2;
1688         }
1689         static
1690         BOOST_UBLAS_INLINE
1691         size_type index_m (size_type index1, size_type /* index2 */) {
1692             return index1;
1693         }
1694         static
1695         BOOST_UBLAS_INLINE
1696         size_type size_M (size_type /* size_i */, size_type size_j) {
1697             return size_j;
1698         }
1699         static
1700         BOOST_UBLAS_INLINE
1701         size_type size_m (size_type size_i, size_type /* size_j */) {
1702             return size_i;
1703         }
1704     };
1705
1706
1707     template <class Z>
1708     struct basic_full {
1709         typedef Z size_type;
1710
1711         template<class L>
1712         static
1713         BOOST_UBLAS_INLINE
1714         size_type packed_size (L, size_type size_i, size_type size_j) {
1715             return L::storage_size (size_i, size_j);
1716         }
1717
1718         static
1719         BOOST_UBLAS_INLINE
1720         bool zero (size_type /* i */, size_type /* j */) {
1721             return false;
1722         }
1723         static
1724         BOOST_UBLAS_INLINE
1725         bool one (size_type /* i */, size_type /* j */) {
1726             return false;
1727         }
1728         static
1729         BOOST_UBLAS_INLINE
1730         bool other (size_type /* i */, size_type /* j */) {
1731             return true;
1732         }
1733         // FIXME: this should not be used at all
1734         static
1735         BOOST_UBLAS_INLINE
1736         size_type restrict1 (size_type i, size_type /* j */) {
1737             return i;
1738         }
1739         static
1740         BOOST_UBLAS_INLINE
1741         size_type restrict2 (size_type /* i */, size_type j) {
1742             return j;
1743         }
1744         static
1745         BOOST_UBLAS_INLINE
1746         size_type mutable_restrict1 (size_type i, size_type /* j */) {
1747             return i;
1748         }
1749         static
1750         BOOST_UBLAS_INLINE
1751         size_type mutable_restrict2 (size_type /* i */, size_type j) {
1752             return j;
1753         }
1754     };
1755
1756     namespace detail {
1757         template < class L >
1758         struct transposed_structure {
1759             typedef typename L::size_type size_type;
1760
1761             template<class LAYOUT>
1762             static
1763             BOOST_UBLAS_INLINE
1764             size_type packed_size (LAYOUT l, size_type size_i, size_type size_j) {
1765                 return L::packed_size(l, size_j, size_i);
1766             }
1767
1768             static
1769             BOOST_UBLAS_INLINE
1770             bool zero (size_type i, size_type j) {
1771                 return L::zero(j, i);
1772             }
1773             static
1774             BOOST_UBLAS_INLINE
1775             bool one (size_type i, size_type j) {
1776                 return L::one(j, i);
1777             }
1778             static
1779             BOOST_UBLAS_INLINE
1780             bool other (size_type i, size_type j) {
1781                 return L::other(j, i);
1782             }
1783             template<class LAYOUT>
1784             static
1785             BOOST_UBLAS_INLINE
1786             size_type element (LAYOUT /* l */, size_type i, size_type size_i, size_type j, size_type size_j) {
1787                 return L::element(typename LAYOUT::transposed_layout(), j, size_j, i, size_i);
1788             }
1789
1790             static
1791             BOOST_UBLAS_INLINE
1792             size_type restrict1 (size_type i, size_type j, size_type size1, size_type size2) {
1793                 return L::restrict2(j, i, size2, size1);
1794             }
1795             static
1796             BOOST_UBLAS_INLINE
1797             size_type restrict2 (size_type i, size_type j, size_type size1, size_type size2) {
1798                 return L::restrict1(j, i, size2, size1);
1799             }
1800             static
1801             BOOST_UBLAS_INLINE
1802             size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type size2) {
1803                 return L::mutable_restrict2(j, i, size2, size1);
1804             }
1805             static
1806             BOOST_UBLAS_INLINE
1807             size_type mutable_restrict2 (size_type i, size_type j, size_type size1, size_type size2) {
1808                 return L::mutable_restrict1(j, i, size2, size1);
1809             }
1810
1811             static
1812             BOOST_UBLAS_INLINE
1813             size_type global_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) {
1814                 return L::global_restrict2(index2, size2, index1, size1);
1815             }
1816             static
1817             BOOST_UBLAS_INLINE
1818             size_type global_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) {
1819                 return L::global_restrict1(index2, size2, index1, size1);
1820             }
1821             static
1822             BOOST_UBLAS_INLINE
1823             size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) {
1824                 return L::global_mutable_restrict2(index2, size2, index1, size1);
1825             }
1826             static
1827             BOOST_UBLAS_INLINE
1828             size_type global_mutable_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) {
1829                 return L::global_mutable_restrict1(index2, size2, index1, size1);
1830             }
1831         };
1832     }
1833
1834     template <class Z>
1835     struct basic_lower {
1836         typedef Z size_type;
1837         typedef lower_tag triangular_type;
1838
1839         template<class L>
1840         static
1841         BOOST_UBLAS_INLINE
1842         size_type packed_size (L, size_type size_i, size_type size_j) {
1843             return L::triangular_size (size_i, size_j);
1844         }
1845
1846         static
1847         BOOST_UBLAS_INLINE
1848         bool zero (size_type i, size_type j) {
1849             return j > i;
1850         }
1851         static
1852         BOOST_UBLAS_INLINE
1853         bool one (size_type /* i */, size_type /* j */) {
1854             return false;
1855         }
1856         static
1857         BOOST_UBLAS_INLINE
1858         bool other (size_type i, size_type j) {
1859             return j <= i;
1860         }
1861         template<class L>
1862         static
1863         BOOST_UBLAS_INLINE
1864         size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) {
1865             return L::lower_element (i, size_i, j, size_j);
1866         }
1867
1868         // return nearest valid index in column j
1869         static
1870         BOOST_UBLAS_INLINE
1871         size_type restrict1 (size_type i, size_type j, size_type size1, size_type /* size2 */) {
1872             return (std::max)(j, (std::min) (size1, i));
1873         }
1874         // return nearest valid index in row i
1875         static
1876         BOOST_UBLAS_INLINE
1877         size_type restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) {
1878             return (std::max)(size_type(0), (std::min) (i+1, j));
1879         }
1880         // return nearest valid mutable index in column j
1881         static
1882         BOOST_UBLAS_INLINE
1883         size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type /* size2 */) {
1884             return (std::max)(j, (std::min) (size1, i));
1885         }
1886         // return nearest valid mutable index in row i
1887         static
1888         BOOST_UBLAS_INLINE
1889         size_type mutable_restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) {
1890             return (std::max)(size_type(0), (std::min) (i+1, j));
1891         }
1892
1893         // return an index between the first and (1+last) filled row
1894         static
1895         BOOST_UBLAS_INLINE
1896         size_type global_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) {
1897             return (std::max)(size_type(0), (std::min)(size1, index1) );
1898         }
1899         // return an index between the first and (1+last) filled column
1900         static
1901         BOOST_UBLAS_INLINE
1902         size_type global_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) {
1903             return (std::max)(size_type(0), (std::min)(size2, index2) );
1904         }
1905
1906         // return an index between the first and (1+last) filled mutable row
1907         static
1908         BOOST_UBLAS_INLINE
1909         size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) {
1910             return (std::max)(size_type(0), (std::min)(size1, index1) );
1911         }
1912         // return an index between the first and (1+last) filled mutable column
1913         static
1914         BOOST_UBLAS_INLINE
1915         size_type global_mutable_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) {
1916             return (std::max)(size_type(0), (std::min)(size2, index2) );
1917         }
1918     };
1919
1920     // the first row only contains a single 1. Thus it is not stored.
1921     template <class Z>
1922     struct basic_unit_lower : public basic_lower<Z> {
1923         typedef Z size_type;
1924         typedef unit_lower_tag triangular_type;
1925
1926         template<class L>
1927         static
1928         BOOST_UBLAS_INLINE
1929         size_type packed_size (L, size_type size_i, size_type size_j) {
1930             // Zero size strict triangles are bad at this point
1931             BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0, bad_index ());
1932             return L::triangular_size (size_i - 1, size_j - 1);
1933         }
1934
1935         static
1936         BOOST_UBLAS_INLINE
1937         bool one (size_type i, size_type j) {
1938             return j == i;
1939         }
1940         static
1941         BOOST_UBLAS_INLINE
1942         bool other (size_type i, size_type j) {
1943             return j < i;
1944         }
1945         template<class L>
1946         static
1947         BOOST_UBLAS_INLINE
1948         size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) {
1949             // Zero size strict triangles are bad at this point
1950             BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0 && i != 0, bad_index ());
1951             return L::lower_element (i-1, size_i - 1, j, size_j - 1);
1952         }
1953
1954         static
1955         BOOST_UBLAS_INLINE
1956         size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type /* size2 */) {
1957             return (std::max)(j+1, (std::min) (size1, i));
1958         }
1959         static
1960         BOOST_UBLAS_INLINE
1961         size_type mutable_restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) {
1962             return (std::max)(size_type(0), (std::min) (i, j));
1963         }
1964
1965         // return an index between the first and (1+last) filled mutable row
1966         static
1967         BOOST_UBLAS_INLINE
1968         size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) {
1969             return (std::max)(size_type(1), (std::min)(size1, index1) );
1970         }
1971         // return an index between the first and (1+last) filled mutable column
1972         static
1973         BOOST_UBLAS_INLINE
1974         size_type global_mutable_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) {
1975             BOOST_UBLAS_CHECK( size2 >= 1 , external_logic() );
1976             return (std::max)(size_type(0), (std::min)(size2-1, index2) );
1977         }
1978     };
1979
1980     // the first row only contains no element. Thus it is not stored.
1981     template <class Z>
1982     struct basic_strict_lower : public basic_unit_lower<Z> {
1983         typedef Z size_type;
1984         typedef strict_lower_tag triangular_type;
1985
1986         template<class L>
1987         static
1988         BOOST_UBLAS_INLINE
1989         size_type packed_size (L, size_type size_i, size_type size_j) {
1990             // Zero size strict triangles are bad at this point
1991             BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0, bad_index ());
1992             return L::triangular_size (size_i - 1, size_j - 1);
1993         }
1994
1995         static
1996         BOOST_UBLAS_INLINE
1997         bool zero (size_type i, size_type j) {
1998             return j >= i;
1999         }
2000         static
2001         BOOST_UBLAS_INLINE
2002         bool one (size_type /*i*/, size_type /*j*/) {
2003             return false;
2004         }
2005         static
2006         BOOST_UBLAS_INLINE
2007         bool other (size_type i, size_type j) {
2008             return j < i;
2009         }
2010         template<class L>
2011         static
2012         BOOST_UBLAS_INLINE
2013         size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) {
2014             // Zero size strict triangles are bad at this point
2015             BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0 && i != 0, bad_index ());
2016             return L::lower_element (i-1, size_i - 1, j, size_j - 1);
2017         }
2018
2019         static
2020         BOOST_UBLAS_INLINE
2021         size_type restrict1 (size_type i, size_type j, size_type size1, size_type size2) {
2022             return basic_unit_lower<Z>::mutable_restrict1(i, j, size1, size2);
2023         }
2024         static
2025         BOOST_UBLAS_INLINE
2026         size_type restrict2 (size_type i, size_type j, size_type size1, size_type size2) {
2027             return basic_unit_lower<Z>::mutable_restrict2(i, j, size1, size2);
2028         }
2029
2030         // return an index between the first and (1+last) filled row
2031         static
2032         BOOST_UBLAS_INLINE
2033         size_type global_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) {
2034             return basic_unit_lower<Z>::global_mutable_restrict1(index1, size1, index2, size2);
2035         }
2036         // return an index between the first and (1+last) filled column
2037         static
2038         BOOST_UBLAS_INLINE
2039         size_type global_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) {
2040             return basic_unit_lower<Z>::global_mutable_restrict2(index1, size1, index2, size2);
2041         }
2042     };
2043
2044
2045     template <class Z>
2046     struct basic_upper : public detail::transposed_structure<basic_lower<Z> >
2047     {
2048         typedef upper_tag triangular_type;
2049     };
2050
2051     template <class Z>
2052     struct basic_unit_upper : public detail::transposed_structure<basic_unit_lower<Z> >
2053     {
2054         typedef unit_upper_tag triangular_type;
2055     };
2056
2057     template <class Z>
2058     struct basic_strict_upper : public detail::transposed_structure<basic_strict_lower<Z> >
2059     {
2060         typedef strict_upper_tag triangular_type;
2061     };
2062
2063
2064 }}}
2065
2066 #endif