boost/numeric/ublas/vector_expression.hpp

   1 //
   2 //  Copyright (c) 2000-2002
   3 //  Joerg Walter, Mathias Koch
   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_VECTOR_EXPRESSION_
  14 #define _BOOST_UBLAS_VECTOR_EXPRESSION_
  15
  16 #include <boost/numeric/ublas/expression_types.hpp>
  17
  18
  19 // Expression templates based on ideas of Todd Veldhuizen and Geoffrey Furnish
  20 // Iterators based on ideas of Jeremy Siek
  21 //
  22 // Classes that model the Vector Expression concept
  23
  24 namespace boost { namespace numeric { namespace ublas {
  25
  26     template<class E>
  27     class vector_reference:
  28         public vector_expression<vector_reference<E> > {
  29
  30         typedef vector_reference<E> self_type;
  31     public:
  32 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
  33         using vector_expression<vector_reference<E> >::operator ();
  34 #endif
  35         typedef typename E::size_type size_type;
  36         typedef typename E::difference_type difference_type;
  37         typedef typename E::value_type value_type;
  38         typedef typename E::const_reference const_reference;
  39         typedef typename boost::mpl::if_<boost::is_const<E>,
  40                                           typename E::const_reference,
  41                                           typename E::reference>::type reference;
  42         typedef E referred_type;
  43         typedef const self_type const_closure_type;
  44         typedef self_type closure_type;
  45         typedef typename E::storage_category storage_category;
  46
  47         // Construction and destruction
  48         BOOST_UBLAS_INLINE
  49         explicit vector_reference (referred_type &e):
  50             e_ (e) {}
  51
  52         // Accessors
  53         BOOST_UBLAS_INLINE
  54         size_type size () const {
  55             return expression ().size ();
  56         }
  57
  58     public:
  59         // Expression accessors - const correct
  60         BOOST_UBLAS_INLINE
  61         const referred_type &expression () const {
  62             return e_;
  63         }
  64         BOOST_UBLAS_INLINE
  65         referred_type &expression () {
  66             return e_;
  67         }
  68
  69     public:
  70         // Element access
  71 #ifndef BOOST_UBLAS_REFERENCE_CONST_MEMBER
  72         BOOST_UBLAS_INLINE
  73         const_reference operator () (size_type i) const {
  74             return expression () (i);
  75         }
  76         BOOST_UBLAS_INLINE
  77         reference operator () (size_type i) {
  78             return expression () (i);
  79         }
  80
  81         BOOST_UBLAS_INLINE
  82         const_reference operator [] (size_type i) const {
  83             return expression () [i];
  84         }
  85         BOOST_UBLAS_INLINE
  86         reference operator [] (size_type i) {
  87             return expression () [i];
  88         }
  89 #else
  90         BOOST_UBLAS_INLINE
  91         reference operator () (size_type i) const {
  92             return expression () (i);
  93         }
  94
  95         BOOST_UBLAS_INLINE
  96         reference operator [] (size_type i) const {
  97             return expression () [i];
  98         }
  99 #endif
 100
 101         // Assignment
 102         BOOST_UBLAS_INLINE
 103         vector_reference &operator = (const vector_reference &v) {
 104             expression ().operator = (v);
 105             return *this;
 106         }
 107         template<class AE>
 108         BOOST_UBLAS_INLINE
 109         vector_reference &operator = (const vector_expression<AE> &ae) {
 110             expression ().operator = (ae);
 111             return *this;
 112         }
 113         template<class AE>
 114         BOOST_UBLAS_INLINE
 115         vector_reference &assign (const vector_expression<AE> &ae) {
 116             expression ().assign (ae);
 117             return *this;
 118         }
 119         template<class AE>
 120         BOOST_UBLAS_INLINE
 121         vector_reference &operator += (const vector_expression<AE> &ae) {
 122             expression ().operator += (ae);
 123             return *this;
 124         }
 125         template<class AE>
 126         BOOST_UBLAS_INLINE
 127         vector_reference &plus_assign (const vector_expression<AE> &ae) {
 128             expression ().plus_assign (ae);
 129             return *this;
 130         }
 131         template<class AE>
 132         BOOST_UBLAS_INLINE
 133         vector_reference &operator -= (const vector_expression<AE> &ae) {
 134             expression ().operator -= (ae);
 135             return *this;
 136         }
 137         template<class AE>
 138         BOOST_UBLAS_INLINE
 139         vector_reference &minus_assign (const vector_expression<AE> &ae) {
 140             expression ().minus_assign (ae);
 141             return *this;
 142         }
 143         template<class AT>
 144         BOOST_UBLAS_INLINE
 145         vector_reference &operator *= (const AT &at) {
 146             expression ().operator *= (at);
 147             return *this;
 148         }
 149         template<class AT>
 150         BOOST_UBLAS_INLINE
 151         vector_reference &operator /= (const AT &at) {
 152             expression ().operator /= (at);
 153             return *this;
 154         }
 155
 156         // Swapping
 157         BOOST_UBLAS_INLINE
 158         void swap (vector_reference &v) {
 159             expression ().swap (v.expression ());
 160         }
 161
 162         // Closure comparison
 163         BOOST_UBLAS_INLINE
 164         bool same_closure (const vector_reference &vr) const {
 165             return &(*this).e_ == &vr.e_;
 166         }
 167
 168         // Iterator types
 169         typedef typename E::const_iterator const_iterator;
 170         typedef typename boost::mpl::if_<boost::is_const<E>,
 171                                           typename E::const_iterator,
 172                                           typename E::iterator>::type iterator;
 173
 174         // Element lookup
 175         BOOST_UBLAS_INLINE
 176         const_iterator find (size_type i) const {
 177             return expression ().find (i);
 178         }
 179         BOOST_UBLAS_INLINE
 180         iterator find (size_type i) {
 181             return expression ().find (i);
 182         }
 183
 184         // Iterator is the iterator of the referenced expression.
 185
 186         BOOST_UBLAS_INLINE
 187         const_iterator begin () const {
 188             return expression ().begin ();
 189         }
 190         BOOST_UBLAS_INLINE
 191         const_iterator end () const {
 192             return expression ().end ();
 193         }
 194
 195         BOOST_UBLAS_INLINE
 196         iterator begin () {
 197             return expression ().begin ();
 198         }
 199         BOOST_UBLAS_INLINE
 200         iterator end () {
 201             return expression ().end ();
 202         }
 203
 204         // Reverse iterator
 205         typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
 206         typedef reverse_iterator_base<iterator> reverse_iterator;
 207
 208         BOOST_UBLAS_INLINE
 209         const_reverse_iterator rbegin () const {
 210             return const_reverse_iterator (end ());
 211         }
 212         BOOST_UBLAS_INLINE
 213         const_reverse_iterator rend () const {
 214             return const_reverse_iterator (begin ());
 215         }
 216         BOOST_UBLAS_INLINE
 217         reverse_iterator rbegin () {
 218             return reverse_iterator (end ());
 219         }
 220         BOOST_UBLAS_INLINE
 221         reverse_iterator rend () {
 222             return reverse_iterator (begin ());
 223         }
 224
 225     private:
 226         referred_type &e_;
 227     };
 228
 229
 230     template<class E, class F>
 231     class vector_unary:
 232         public vector_expression<vector_unary<E, F> > {
 233
 234         typedef F functor_type;
 235         typedef typename boost::mpl::if_<boost::is_same<F, scalar_identity<typename E::value_type> >,
 236                                           E,
 237                                           const E>::type expression_type;
 238         typedef typename boost::mpl::if_<boost::is_const<expression_type>,
 239                                           typename E::const_closure_type,
 240                                           typename E::closure_type>::type expression_closure_type;
 241         typedef vector_unary<E, F> self_type;
 242     public:
 243 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
 244         using vector_expression<vector_unary<E, F> >::operator ();
 245 #endif
 246         typedef typename E::size_type size_type;
 247         typedef typename E::difference_type difference_type;
 248         typedef typename F::result_type value_type;
 249         typedef value_type const_reference;
 250         typedef typename boost::mpl::if_<boost::is_same<F, scalar_identity<value_type> >,
 251                                           typename E::reference,
 252                                           value_type>::type reference;
 253         typedef const self_type const_closure_type;
 254         typedef self_type closure_type;
 255         typedef unknown_storage_tag storage_category;
 256
 257         // Construction and destruction
 258         BOOST_UBLAS_INLINE
 259         // May be used as mutable expression.
 260         explicit vector_unary (expression_type &e):
 261             e_ (e) {}
 262
 263         // Accessors
 264         BOOST_UBLAS_INLINE
 265         size_type size () const {
 266             return e_.size ();
 267         }
 268
 269     public:
 270         // Expression accessors
 271         BOOST_UBLAS_INLINE
 272         const expression_closure_type &expression () const {
 273             return e_;
 274         }
 275
 276     public:
 277         // Element access
 278         BOOST_UBLAS_INLINE
 279         const_reference operator () (size_type i) const {
 280             return functor_type::apply (e_ (i));
 281         }
 282         BOOST_UBLAS_INLINE
 283         reference operator () (size_type i) {
 284             BOOST_STATIC_ASSERT ((boost::is_same<functor_type, scalar_identity<value_type > >::value));
 285             return e_ (i);
 286         }
 287
 288         BOOST_UBLAS_INLINE
 289         const_reference operator [] (size_type i) const {
 290             return functor_type::apply (e_ [i]);
 291         }
 292         BOOST_UBLAS_INLINE
 293         reference operator [] (size_type i) {
 294             BOOST_STATIC_ASSERT ((boost::is_same<functor_type, scalar_identity<value_type > >::value));
 295             return e_ [i];
 296         }
 297
 298         // Closure comparison
 299         BOOST_UBLAS_INLINE
 300         bool same_closure (const vector_unary &vu) const {
 301             return (*this).expression ().same_closure (vu.expression ());
 302         }
 303
 304         // Iterator types
 305     private:
 306         typedef typename E::const_iterator const_subiterator_type;
 307         typedef const value_type *const_pointer;
 308
 309     public:
 310 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
 311         typedef indexed_const_iterator<const_closure_type, typename const_subiterator_type::iterator_category> const_iterator;
 312         typedef const_iterator iterator;
 313 #else
 314         class const_iterator;
 315         typedef const_iterator iterator;
 316 #endif
 317
 318         // Element lookup
 319         BOOST_UBLAS_INLINE
 320         const_iterator find (size_type i) const {
 321 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
 322             const_subiterator_type it (e_.find (i));
 323             return const_iterator (*this, it.index ());
 324 #else
 325             return const_iterator (*this, e_.find (i));
 326 #endif
 327         }
 328
 329         // Iterator enhances the iterator of the referenced expression
 330         // with the unary functor.
 331
 332 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
 333         class const_iterator:
 334             public container_const_reference<vector_unary>,
 335             public iterator_base_traits<typename E::const_iterator::iterator_category>::template
 336                         iterator_base<const_iterator, value_type>::type {
 337         public:
 338             typedef typename E::const_iterator::iterator_category iterator_category;
 339             typedef typename vector_unary::difference_type difference_type;
 340             typedef typename vector_unary::value_type value_type;
 341             typedef typename vector_unary::const_reference reference;
 342             typedef typename vector_unary::const_pointer pointer;
 343
 344             // Construction and destruction
 345             BOOST_UBLAS_INLINE
 346             const_iterator ():
 347                 container_const_reference<self_type> (), it_ () {}
 348             BOOST_UBLAS_INLINE
 349             const_iterator (const self_type &vu, const const_subiterator_type &it):
 350                 container_const_reference<self_type> (vu), it_ (it) {}
 351
 352             // Arithmetic
 353             BOOST_UBLAS_INLINE
 354             const_iterator &operator ++ () {
 355                 ++ it_;
 356                 return *this;
 357             }
 358             BOOST_UBLAS_INLINE
 359             const_iterator &operator -- () {
 360                 -- it_;
 361                 return *this;
 362             }
 363             BOOST_UBLAS_INLINE
 364             const_iterator &operator += (difference_type n) {
 365                 it_ += n;
 366                 return *this;
 367             }
 368             BOOST_UBLAS_INLINE
 369             const_iterator &operator -= (difference_type n) {
 370                 it_ -= n;
 371                 return *this;
 372             }
 373             BOOST_UBLAS_INLINE
 374             difference_type operator - (const const_iterator &it) const {
 375                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
 376                 return it_ - it.it_;
 377             }
 378
 379             // Dereference
 380             BOOST_UBLAS_INLINE
 381             const_reference operator * () const {
 382                 return functor_type::apply (*it_);
 383             }
 384             BOOST_UBLAS_INLINE
 385             const_reference operator [] (difference_type n) const {
 386                 return *(*this + n);
 387             }
 388
 389             // Index
 390             BOOST_UBLAS_INLINE
 391             size_type index () const {
 392                 return it_.index ();
 393             }
 394
 395             // Assignment
 396             BOOST_UBLAS_INLINE
 397             const_iterator &operator = (const const_iterator &it) {
 398                 container_const_reference<self_type>::assign (&it ());
 399                 it_ = it.it_;
 400                 return *this;
 401             }
 402
 403             // Comparison
 404             BOOST_UBLAS_INLINE
 405             bool operator == (const const_iterator &it) const {
 406                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
 407                 return it_ == it.it_;
 408             }
 409             BOOST_UBLAS_INLINE
 410             bool operator < (const const_iterator &it) const {
 411                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
 412                 return it_ < it.it_;
 413             }
 414
 415         private:
 416             const_subiterator_type it_;
 417         };
 418 #endif
 419
 420         BOOST_UBLAS_INLINE
 421         const_iterator begin () const {
 422             return find (0);
 423         }
 424         BOOST_UBLAS_INLINE
 425         const_iterator end () const {
 426             return find (size ());
 427         }
 428
 429         // Reverse iterator
 430         typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
 431
 432         BOOST_UBLAS_INLINE
 433         const_reverse_iterator rbegin () const {
 434             return const_reverse_iterator (end ());
 435         }
 436         BOOST_UBLAS_INLINE
 437         const_reverse_iterator rend () const {
 438             return const_reverse_iterator (begin ());
 439         }
 440
 441     private:
 442         expression_closure_type e_;
 443     };
 444
 445     template<class E, class F>
 446     struct vector_unary_traits {
 447         typedef vector_unary<E, F> expression_type;
 448 //FIXME
 449 // #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
 450         typedef expression_type result_type;
 451 // #else
 452 //         typedef typename E::vector_temporary_type result_type;
 453 // #endif
 454     };
 455
 456     // (- v) [i] = - v [i]
 457     template<class E>
 458     BOOST_UBLAS_INLINE
 459     typename vector_unary_traits<E, scalar_negate<typename E::value_type> >::result_type
 460     operator - (const vector_expression<E> &e) {
 461         typedef typename vector_unary_traits<E, scalar_negate<typename E::value_type> >::expression_type expression_type;
 462         return expression_type (e ());
 463     }
 464
 465     // (conj v) [i] = conj (v [i])
 466     template<class E>
 467     BOOST_UBLAS_INLINE
 468     typename vector_unary_traits<E, scalar_conj<typename E::value_type> >::result_type
 469     conj (const vector_expression<E> &e) {
 470         typedef typename vector_unary_traits<E, scalar_conj<typename E::value_type> >::expression_type expression_type;
 471         return expression_type (e ());
 472     }
 473
 474     // (real v) [i] = real (v [i])
 475     template<class E>
 476     BOOST_UBLAS_INLINE
 477     typename vector_unary_traits<E, scalar_real<typename E::value_type> >::result_type
 478     real (const vector_expression<E> &e) {
 479         typedef typename vector_unary_traits<E, scalar_real<typename E::value_type> >::expression_type expression_type;
 480         return expression_type (e ());
 481     }
 482
 483     // (imag v) [i] = imag (v [i])
 484     template<class E>
 485     BOOST_UBLAS_INLINE
 486     typename vector_unary_traits<E, scalar_imag<typename E::value_type> >::result_type
 487     imag (const vector_expression<E> &e) {
 488         typedef typename vector_unary_traits<E, scalar_imag<typename E::value_type> >::expression_type expression_type;
 489         return expression_type (e ());
 490     }
 491
 492     // (trans v) [i] = v [i]
 493     template<class E>
 494     BOOST_UBLAS_INLINE
 495     typename vector_unary_traits<const E, scalar_identity<typename E::value_type> >::result_type
 496     trans (const vector_expression<E> &e) {
 497         typedef typename vector_unary_traits<const E, scalar_identity<typename E::value_type> >::expression_type expression_type;
 498         return expression_type (e ());
 499     }
 500     template<class E>
 501     BOOST_UBLAS_INLINE
 502     typename vector_unary_traits<E, scalar_identity<typename E::value_type> >::result_type
 503     trans (vector_expression<E> &e) {
 504         typedef typename vector_unary_traits<E, scalar_identity<typename E::value_type> >::expression_type expression_type;
 505         return expression_type (e ());
 506     }
 507
 508     // (herm v) [i] = conj (v [i])
 509     template<class E>
 510     BOOST_UBLAS_INLINE
 511     typename vector_unary_traits<E, scalar_conj<typename E::value_type> >::result_type
 512     herm (const vector_expression<E> &e) {
 513         typedef typename vector_unary_traits<E, scalar_conj<typename E::value_type> >::expression_type expression_type;
 514         return expression_type (e ());
 515     }
 516
 517     template<class E1, class E2, class F>
 518     class vector_binary:
 519         public vector_expression<vector_binary<E1, E2, F> > {
 520
 521         typedef E1 expression1_type;
 522         typedef E2 expression2_type;
 523         typedef F functor_type;
 524         typedef typename E1::const_closure_type expression1_closure_type;
 525         typedef typename E2::const_closure_type expression2_closure_type;
 526         typedef vector_binary<E1, E2, F> self_type;
 527     public:
 528 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
 529         using vector_expression<vector_binary<E1, E2, F> >::operator ();
 530 #endif
 531         typedef typename promote_traits<typename E1::size_type, typename E2::size_type>::promote_type size_type;
 532         typedef typename promote_traits<typename E1::difference_type, typename E2::difference_type>::promote_type difference_type;
 533         typedef typename F::result_type value_type;
 534         typedef value_type const_reference;
 535         typedef const_reference reference;
 536         typedef const self_type const_closure_type;
 537         typedef const_closure_type closure_type;
 538         typedef unknown_storage_tag storage_category;
 539
 540         // Construction and destruction
 541         BOOST_UBLAS_INLINE
 542         vector_binary (const expression1_type &e1, const expression2_type &e2):
 543             e1_ (e1), e2_ (e2) {}
 544
 545         // Accessors
 546         BOOST_UBLAS_INLINE
 547         size_type size () const {
 548             return BOOST_UBLAS_SAME (e1_.size (), e2_.size ());
 549         }
 550
 551     private:
 552         // Accessors
 553         BOOST_UBLAS_INLINE
 554         const expression1_closure_type &expression1 () const {
 555             return e1_;
 556         }
 557         BOOST_UBLAS_INLINE
 558         const expression2_closure_type &expression2 () const {
 559             return e2_;
 560         }
 561
 562     public:
 563         // Element access
 564         BOOST_UBLAS_INLINE
 565         const_reference operator () (size_type i) const {
 566             return functor_type::apply (e1_ (i), e2_ (i));
 567         }
 568
 569         BOOST_UBLAS_INLINE
 570         const_reference operator [] (size_type i) const {
 571             return functor_type::apply (e1_ [i], e2_ [i]);
 572         }
 573
 574         // Closure comparison
 575         BOOST_UBLAS_INLINE
 576         bool same_closure (const vector_binary &vb) const {
 577             return (*this).expression1 ().same_closure (vb.expression1 ()) &&
 578                    (*this).expression2 ().same_closure (vb.expression2 ());
 579         }
 580
 581         // Iterator types
 582     private:
 583         typedef typename E1::const_iterator const_subiterator1_type;
 584         typedef typename E2::const_iterator const_subiterator2_type;
 585         typedef const value_type *const_pointer;
 586
 587     public:
 588 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
 589         typedef typename iterator_restrict_traits<typename const_subiterator1_type::iterator_category,
 590                                                   typename const_subiterator2_type::iterator_category>::iterator_category iterator_category;
 591         typedef indexed_const_iterator<const_closure_type, iterator_category> const_iterator;
 592         typedef const_iterator iterator;
 593 #else
 594         class const_iterator;
 595         typedef const_iterator iterator;
 596 #endif
 597
 598         // Element lookup
 599         BOOST_UBLAS_INLINE
 600         const_iterator find (size_type i) const {
 601             const_subiterator1_type it1 (e1_.find (i));
 602             const_subiterator1_type it1_end (e1_.find (size ()));
 603             const_subiterator2_type it2 (e2_.find (i));
 604             const_subiterator2_type it2_end (e2_.find (size ()));
 605             i = (std::min) (it1 != it1_end ? it1.index () : size (),
 606                           it2 != it2_end ? it2.index () : size ());
 607 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
 608             return const_iterator (*this, i);
 609 #else
 610             return const_iterator (*this, i, it1, it1_end, it2, it2_end);
 611 #endif
 612         }
 613
 614         // Iterator merges the iterators of the referenced expressions and
 615         // enhances them with the binary functor.
 616
 617 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
 618         class const_iterator:
 619             public container_const_reference<vector_binary>,
 620             public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
 621                                                                           typename E2::const_iterator::iterator_category>::iterator_category>::template
 622                         iterator_base<const_iterator, value_type>::type {
 623         public:
 624             typedef typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
 625                                                       typename E2::const_iterator::iterator_category>::iterator_category iterator_category;
 626             typedef typename vector_binary::difference_type difference_type;
 627             typedef typename vector_binary::value_type value_type;
 628             typedef typename vector_binary::const_reference reference;
 629             typedef typename vector_binary::const_pointer pointer;
 630
 631             // Construction and destruction
 632             BOOST_UBLAS_INLINE
 633             const_iterator ():
 634                 container_const_reference<self_type> (), i_ (), it1_ (), it1_end_ (), it2_ (), it2_end_ () {}
 635             BOOST_UBLAS_INLINE
 636             const_iterator (const self_type &vb, size_type i,
 637                             const const_subiterator1_type &it1, const const_subiterator1_type &it1_end,
 638                             const const_subiterator2_type &it2, const const_subiterator2_type &it2_end):
 639                 container_const_reference<self_type> (vb), i_ (i), it1_ (it1), it1_end_ (it1_end), it2_ (it2), it2_end_ (it2_end) {}
 640
 641         private:
 642             // Dense specializations
 643             BOOST_UBLAS_INLINE
 644             void increment (dense_random_access_iterator_tag) {
 645                 ++ i_; ++ it1_; ++ it2_;
 646             }
 647             BOOST_UBLAS_INLINE
 648             void decrement (dense_random_access_iterator_tag) {
 649                 -- i_; -- it1_; -- it2_;
 650             }
 651             BOOST_UBLAS_INLINE
 652             void increment (dense_random_access_iterator_tag, difference_type n) {
 653                 i_ += n; it1_ += n; it2_ += n;
 654             }
 655             BOOST_UBLAS_INLINE
 656             void decrement (dense_random_access_iterator_tag, difference_type n) {
 657                 i_ -= n; it1_ -= n; it2_ -= n;
 658             }
 659             BOOST_UBLAS_INLINE
 660             value_type dereference (dense_random_access_iterator_tag) const {
 661                 return functor_type::apply (*it1_, *it2_);
 662             }
 663
 664             // Packed specializations
 665             BOOST_UBLAS_INLINE
 666             void increment (packed_random_access_iterator_tag) {
 667                 if (it1_ != it1_end_)
 668                     if (it1_.index () <= i_)
 669                         ++ it1_;
 670                 if (it2_ != it2_end_)
 671                     if (it2_.index () <= i_)
 672                         ++ it2_;
 673                 ++ i_;
 674             }
 675             BOOST_UBLAS_INLINE
 676             void decrement (packed_random_access_iterator_tag) {
 677                 if (it1_ != it1_end_)
 678                     if (i_ <= it1_.index ())
 679                         -- it1_;
 680                 if (it2_ != it2_end_)
 681                     if (i_ <= it2_.index ())
 682                         -- it2_;
 683                 -- i_;
 684             }
 685             BOOST_UBLAS_INLINE
 686             void increment (packed_random_access_iterator_tag, difference_type n) {
 687                 while (n > 0) {
 688                     increment (packed_random_access_iterator_tag ());
 689                     --n;
 690                 }
 691                 while (n < 0) {
 692                     decrement (packed_random_access_iterator_tag ());
 693                     ++n;
 694                 }
 695             }
 696             BOOST_UBLAS_INLINE
 697             void decrement (packed_random_access_iterator_tag, difference_type n) {
 698                 while (n > 0) {
 699                     decrement (packed_random_access_iterator_tag ());
 700                     --n;
 701                 }
 702                 while (n < 0) {
 703                     increment (packed_random_access_iterator_tag ());
 704                     ++n;
 705                 }
 706             }
 707             BOOST_UBLAS_INLINE
 708             value_type dereference (packed_random_access_iterator_tag) const {
 709                 value_type t1 = value_type/*zero*/();
 710                 if (it1_ != it1_end_)
 711                     if (it1_.index () == i_)
 712                         t1 = *it1_;
 713                 value_type t2 = value_type/*zero*/();
 714                 if (it2_ != it2_end_)
 715                     if (it2_.index () == i_)
 716                         t2 = *it2_;
 717                 return functor_type::apply (t1, t2);
 718             }
 719
 720             // Sparse specializations
 721             BOOST_UBLAS_INLINE
 722             void increment (sparse_bidirectional_iterator_tag) {
 723                 size_type index1 = (*this) ().size ();
 724                 if (it1_ != it1_end_) {
 725                     if  (it1_.index () <= i_)
 726                         ++ it1_;
 727                     if (it1_ != it1_end_)
 728                         index1 = it1_.index ();
 729                 }
 730                 size_type index2 = (*this) ().size ();
 731                 if (it2_ != it2_end_) {
 732                     if (it2_.index () <= i_)
 733                         ++ it2_;
 734                     if (it2_ != it2_end_)
 735                         index2 = it2_.index ();
 736                 }
 737                 i_ = (std::min) (index1, index2);
 738             }
 739             BOOST_UBLAS_INLINE
 740             void decrement (sparse_bidirectional_iterator_tag) {
 741                 size_type index1 = (*this) ().size ();
 742                 if (it1_ != it1_end_) {
 743                     if (i_ <= it1_.index ())
 744                         -- it1_;
 745                     if (it1_ != it1_end_)
 746                         index1 = it1_.index ();
 747                 }
 748                 size_type index2 = (*this) ().size ();
 749                 if (it2_ != it2_end_) {
 750                     if (i_ <= it2_.index ())
 751                         -- it2_;
 752                     if (it2_ != it2_end_)
 753                         index2 = it2_.index ();
 754                 }
 755                 i_ = (std::max) (index1, index2);
 756             }
 757             BOOST_UBLAS_INLINE
 758             void increment (sparse_bidirectional_iterator_tag, difference_type n) {
 759                 while (n > 0) {
 760                     increment (sparse_bidirectional_iterator_tag ());
 761                     --n;
 762                 }
 763                 while (n < 0) {
 764                     decrement (sparse_bidirectional_iterator_tag ());
 765                     ++n;
 766                 }
 767             }
 768             BOOST_UBLAS_INLINE
 769             void decrement (sparse_bidirectional_iterator_tag, difference_type n) {
 770                 while (n > 0) {
 771                     decrement (sparse_bidirectional_iterator_tag ());
 772                     --n;
 773                 }
 774                 while (n < 0) {
 775                     increment (sparse_bidirectional_iterator_tag ());
 776                     ++n;
 777                 }
 778             }
 779             BOOST_UBLAS_INLINE
 780             value_type dereference (sparse_bidirectional_iterator_tag) const {
 781                 value_type t1 = value_type/*zero*/();
 782                 if (it1_ != it1_end_)
 783                     if (it1_.index () == i_)
 784                         t1 = *it1_;
 785                 value_type t2 = value_type/*zero*/();
 786                 if (it2_ != it2_end_)
 787                     if (it2_.index () == i_)
 788                         t2 = *it2_;
 789                 return functor_type::apply (t1, t2);
 790             }
 791
 792         public:
 793             // Arithmetic
 794             BOOST_UBLAS_INLINE
 795             const_iterator &operator ++ () {
 796                 increment (iterator_category ());
 797                 return *this;
 798             }
 799             BOOST_UBLAS_INLINE
 800             const_iterator &operator -- () {
 801                 decrement (iterator_category ());
 802                 return *this;
 803             }
 804             BOOST_UBLAS_INLINE
 805             const_iterator &operator += (difference_type n) {
 806                 increment (iterator_category (), n);
 807                 return *this;
 808             }
 809             BOOST_UBLAS_INLINE
 810             const_iterator &operator -= (difference_type n) {
 811                 decrement (iterator_category (), n);
 812                 return *this;
 813             }
 814             BOOST_UBLAS_INLINE
 815             difference_type operator - (const const_iterator &it) const {
 816                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
 817                 return index () - it.index ();
 818             }
 819
 820             // Dereference
 821             BOOST_UBLAS_INLINE
 822             const_reference operator * () const {
 823                 return dereference (iterator_category ());
 824             }
 825             BOOST_UBLAS_INLINE
 826             const_reference operator [] (difference_type n) const {
 827                 return *(*this + n);
 828             }
 829
 830             // Index
 831             BOOST_UBLAS_INLINE
 832             size_type index () const {
 833                 return i_;
 834             }
 835
 836             // Assignment
 837             BOOST_UBLAS_INLINE
 838             const_iterator &operator = (const const_iterator &it) {
 839                 container_const_reference<self_type>::assign (&it ());
 840                 i_ = it.i_;
 841                 it1_ = it.it1_;
 842                 it1_end_ = it.it1_end_;
 843                 it2_ = it.it2_;
 844                 it2_end_ = it.it2_end_;
 845                 return *this;
 846             }
 847
 848             // Comparison
 849             BOOST_UBLAS_INLINE
 850             bool operator == (const const_iterator &it) const {
 851                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
 852                 return index () == it.index ();
 853             }
 854             BOOST_UBLAS_INLINE
 855             bool operator < (const const_iterator &it) const {
 856                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
 857                 return index () < it.index ();
 858             }
 859
 860         private:
 861             size_type i_;
 862             const_subiterator1_type it1_;
 863             const_subiterator1_type it1_end_;
 864             const_subiterator2_type it2_;
 865             const_subiterator2_type it2_end_;
 866         };
 867 #endif
 868
 869         BOOST_UBLAS_INLINE
 870         const_iterator begin () const {
 871             return find (0);
 872         }
 873         BOOST_UBLAS_INLINE
 874         const_iterator end () const {
 875             return find (size ());
 876         }
 877
 878         // Reverse iterator
 879         typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
 880
 881         BOOST_UBLAS_INLINE
 882         const_reverse_iterator rbegin () const {
 883             return const_reverse_iterator (end ());
 884         }
 885         BOOST_UBLAS_INLINE
 886         const_reverse_iterator rend () const {
 887             return const_reverse_iterator (begin ());
 888         }
 889
 890     private:
 891         expression1_closure_type e1_;
 892         expression2_closure_type e2_;
 893     };
 894
 895     template<class E1, class E2, class F>
 896     struct vector_binary_traits {
 897         typedef vector_binary<E1, E2, F> expression_type;
 898 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
 899         typedef expression_type result_type;
 900 #else
 901         typedef typename E1::vector_temporary_type result_type;
 902 #endif
 903     };
 904
 905     // (v1 + v2) [i] = v1 [i] + v2 [i]
 906     template<class E1, class E2>
 907     BOOST_UBLAS_INLINE
 908     typename vector_binary_traits<E1, E2, scalar_plus<typename E1::value_type,
 909                                                       typename E2::value_type> >::result_type
 910     operator + (const vector_expression<E1> &e1,
 911                 const vector_expression<E2> &e2) {
 912         typedef typename vector_binary_traits<E1, E2, scalar_plus<typename E1::value_type,
 913                                                                               typename E2::value_type> >::expression_type expression_type;
 914         return expression_type (e1 (), e2 ());
 915     }
 916
 917     // (v1 - v2) [i] = v1 [i] - v2 [i]
 918     template<class E1, class E2>
 919     BOOST_UBLAS_INLINE
 920     typename vector_binary_traits<E1, E2, scalar_minus<typename E1::value_type,
 921                                                        typename E2::value_type> >::result_type
 922     operator - (const vector_expression<E1> &e1,
 923                 const vector_expression<E2> &e2) {
 924         typedef typename vector_binary_traits<E1, E2, scalar_minus<typename E1::value_type,
 925                                                                                typename E2::value_type> >::expression_type expression_type;
 926         return expression_type (e1 (), e2 ());
 927     }
 928
 929     // (v1 * v2) [i] = v1 [i] * v2 [i]
 930     template<class E1, class E2>
 931     BOOST_UBLAS_INLINE
 932     typename vector_binary_traits<E1, E2, scalar_multiplies<typename E1::value_type,
 933                                                             typename E2::value_type> >::result_type
 934     element_prod (const vector_expression<E1> &e1,
 935                   const vector_expression<E2> &e2) {
 936         typedef typename vector_binary_traits<E1, E2, scalar_multiplies<typename E1::value_type,
 937                                                                                     typename E2::value_type> >::expression_type expression_type;
 938         return expression_type (e1 (), e2 ());
 939     }
 940
 941     // (v1 / v2) [i] = v1 [i] / v2 [i]
 942     template<class E1, class E2>
 943     BOOST_UBLAS_INLINE
 944     typename vector_binary_traits<E1, E2, scalar_divides<typename E1::value_type,
 945                                                          typename E2::value_type> >::result_type
 946     element_div (const vector_expression<E1> &e1,
 947                  const vector_expression<E2> &e2) {
 948         typedef typename vector_binary_traits<E1, E2, scalar_divides<typename E1::value_type,
 949                                                                                  typename E2::value_type> >::expression_type expression_type;
 950         return expression_type (e1 (), e2 ());
 951     }
 952
 953
 954     template<class E1, class E2, class F>
 955     class vector_binary_scalar1:
 956         public vector_expression<vector_binary_scalar1<E1, E2, F> > {
 957
 958         typedef F functor_type;
 959         typedef E1 expression1_type;
 960         typedef E2 expression2_type;
 961     public:
 962         typedef const E1& expression1_closure_type;
 963         typedef typename E2::const_closure_type expression2_closure_type;
 964     private:
 965         typedef vector_binary_scalar1<E1, E2, F> self_type;
 966     public:
 967 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
 968         using vector_expression<vector_binary_scalar1<E1, E2, F> >::operator ();
 969 #endif
 970         typedef typename E2::size_type size_type;
 971         typedef typename E2::difference_type difference_type;
 972         typedef typename F::result_type value_type;
 973         typedef value_type const_reference;
 974         typedef const_reference reference;
 975         typedef const self_type const_closure_type;
 976         typedef const_closure_type closure_type;
 977         typedef unknown_storage_tag storage_category;
 978
 979         // Construction and destruction
 980         BOOST_UBLAS_INLINE
 981         vector_binary_scalar1 (const expression1_type &e1, const expression2_type &e2):
 982             e1_ (e1), e2_ (e2) {}
 983
 984         // Accessors
 985         BOOST_UBLAS_INLINE
 986         size_type size () const {
 987             return e2_.size ();
 988         }
 989
 990     public:
 991         // Element access
 992         BOOST_UBLAS_INLINE
 993         const_reference operator () (size_type i) const {
 994             return functor_type::apply (e1_, e2_ (i));
 995         }
 996
 997         BOOST_UBLAS_INLINE
 998         const_reference operator [] (size_type i) const {
 999             return functor_type::apply (e1_, e2_ [i]);
1000         }
1001
1002         // Closure comparison
1003         BOOST_UBLAS_INLINE
1004         bool same_closure (const vector_binary_scalar1 &vbs1) const {
1005             return &e1_ == &(vbs1.e1_) &&
1006                    (*this).e2_.same_closure (vbs1.e2_);
1007         }
1008
1009         // Iterator types
1010     private:
1011         typedef expression1_type const_subiterator1_type;
1012         typedef typename expression2_type::const_iterator const_subiterator2_type;
1013         typedef const value_type *const_pointer;
1014
1015     public:
1016 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1017         typedef indexed_const_iterator<const_closure_type, typename const_subiterator2_type::iterator_category> const_iterator;
1018         typedef const_iterator iterator;
1019 #else
1020         class const_iterator;
1021         typedef const_iterator iterator;
1022 #endif
1023
1024         // Element lookup
1025         BOOST_UBLAS_INLINE
1026         const_iterator find (size_type i) const {
1027 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1028             const_subiterator2_type it (e2_.find (i));
1029             return const_iterator (*this, it.index ());
1030 #else
1031             return const_iterator (*this, const_subiterator1_type (e1_), e2_.find (i));
1032 #endif
1033         }
1034
1035         // Iterator enhances the iterator of the referenced vector expression
1036         // with the binary functor.
1037
1038 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1039         class const_iterator:
1040             public container_const_reference<vector_binary_scalar1>,
1041             public iterator_base_traits<typename E2::const_iterator::iterator_category>::template
1042                         iterator_base<const_iterator, value_type>::type {
1043         public:
1044             typedef typename E2::const_iterator::iterator_category iterator_category;
1045             typedef typename vector_binary_scalar1::difference_type difference_type;
1046             typedef typename vector_binary_scalar1::value_type value_type;
1047             typedef typename vector_binary_scalar1::const_reference reference;
1048             typedef typename vector_binary_scalar1::const_pointer pointer;
1049
1050             // Construction and destruction
1051             BOOST_UBLAS_INLINE
1052             const_iterator ():
1053                 container_const_reference<self_type> (), it1_ (), it2_ () {}
1054             BOOST_UBLAS_INLINE
1055             const_iterator (const self_type &vbs, const const_subiterator1_type &it1, const const_subiterator2_type &it2):
1056                 container_const_reference<self_type> (vbs), it1_ (it1), it2_ (it2) {}
1057
1058             // Arithmetic
1059             BOOST_UBLAS_INLINE
1060             const_iterator &operator ++ () {
1061                 ++ it2_;
1062                 return *this;
1063             }
1064             BOOST_UBLAS_INLINE
1065             const_iterator &operator -- () {
1066                 -- it2_;
1067                 return *this;
1068             }
1069             BOOST_UBLAS_INLINE
1070             const_iterator &operator += (difference_type n) {
1071                 it2_ += n;
1072                 return *this;
1073             }
1074             BOOST_UBLAS_INLINE
1075             const_iterator &operator -= (difference_type n) {
1076                 it2_ -= n;
1077                 return *this;
1078             }
1079             BOOST_UBLAS_INLINE
1080             difference_type operator - (const const_iterator &it) const {
1081                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1082                 // FIXME we shouldn't compare floats
1083                 // BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
1084                 return it2_ - it.it2_;
1085             }
1086
1087             // Dereference
1088             BOOST_UBLAS_INLINE
1089             const_reference operator * () const {
1090                 return functor_type::apply (it1_, *it2_);
1091             }
1092             BOOST_UBLAS_INLINE
1093             const_reference operator [] (difference_type n) const {
1094                 return *(*this + n);
1095             }
1096
1097             // Index
1098             BOOST_UBLAS_INLINE
1099             size_type index () const {
1100                 return it2_.index ();
1101             }
1102
1103             // Assignment
1104             BOOST_UBLAS_INLINE
1105             const_iterator &operator = (const const_iterator &it) {
1106                 container_const_reference<self_type>::assign (&it ());
1107                 it1_ = it.it1_;
1108                 it2_ = it.it2_;
1109                 return *this;
1110             }
1111
1112             // Comparison
1113             BOOST_UBLAS_INLINE
1114             bool operator == (const const_iterator &it) const {
1115                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1116                 // FIXME we shouldn't compare floats
1117                 // BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
1118                 return it2_ == it.it2_;
1119             }
1120             BOOST_UBLAS_INLINE
1121             bool operator < (const const_iterator &it) const {
1122                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1123                 // FIXME we shouldn't compare floats
1124                 // BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
1125                 return it2_ < it.it2_;
1126             }
1127
1128         private:
1129             const_subiterator1_type it1_;
1130             const_subiterator2_type it2_;
1131         };
1132 #endif
1133
1134         BOOST_UBLAS_INLINE
1135         const_iterator begin () const {
1136             return find (0);
1137         }
1138         BOOST_UBLAS_INLINE
1139         const_iterator end () const {
1140             return find (size ());
1141         }
1142
1143         // Reverse iterator
1144         typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
1145
1146         BOOST_UBLAS_INLINE
1147         const_reverse_iterator rbegin () const {
1148             return const_reverse_iterator (end ());
1149         }
1150         BOOST_UBLAS_INLINE
1151         const_reverse_iterator rend () const {
1152             return const_reverse_iterator (begin ());
1153         }
1154
1155     private:
1156         expression1_closure_type e1_;
1157         expression2_closure_type e2_;
1158     };
1159
1160     template<class E1, class E2, class F>
1161     struct vector_binary_scalar1_traits {
1162         typedef vector_binary_scalar1<E1, E2, F> expression_type;   // allow E1 to be builtin type
1163 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
1164         typedef expression_type result_type;
1165 #else
1166         typedef typename E2::vector_temporary_type result_type;
1167 #endif
1168     };
1169
1170     // (t * v) [i] = t * v [i]
1171     template<class T1, class E2>
1172     BOOST_UBLAS_INLINE
1173     typename vector_binary_scalar1_traits<const T1, E2, scalar_multiplies<T1, typename E2::value_type> >::result_type
1174     operator * (const T1 &e1,
1175                 const vector_expression<E2> &e2) {
1176         typedef typename vector_binary_scalar1_traits<const T1, E2, scalar_multiplies<T1, typename E2::value_type> >::expression_type expression_type;
1177         return expression_type (e1, e2 ());
1178     }
1179
1180
1181     template<class E1, class E2, class F>
1182     class vector_binary_scalar2:
1183         public vector_expression<vector_binary_scalar2<E1, E2, F> > {
1184
1185         typedef F functor_type;
1186         typedef E1 expression1_type;
1187         typedef E2 expression2_type;
1188         typedef typename E1::const_closure_type expression1_closure_type;
1189         typedef const E2& expression2_closure_type;
1190         typedef vector_binary_scalar2<E1, E2, F> self_type;
1191     public:
1192 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
1193         using vector_expression<vector_binary_scalar2<E1, E2, F> >::operator ();
1194 #endif
1195         typedef typename E1::size_type size_type;
1196         typedef typename E1::difference_type difference_type;
1197         typedef typename F::result_type value_type;
1198         typedef value_type const_reference;
1199         typedef const_reference reference;
1200         typedef const self_type const_closure_type;
1201         typedef const_closure_type closure_type;
1202         typedef unknown_storage_tag storage_category;
1203
1204         // Construction and destruction
1205         BOOST_UBLAS_INLINE
1206         vector_binary_scalar2 (const expression1_type &e1, const expression2_type &e2):
1207             e1_ (e1), e2_ (e2) {}
1208
1209         // Accessors
1210         BOOST_UBLAS_INLINE
1211         size_type size () const {
1212             return e1_.size ();
1213         }
1214
1215     public:
1216         // Element access
1217         BOOST_UBLAS_INLINE
1218         const_reference operator () (size_type i) const {
1219             return functor_type::apply (e1_ (i), e2_);
1220         }
1221
1222         BOOST_UBLAS_INLINE
1223         const_reference operator [] (size_type i) const {
1224             return functor_type::apply (e1_ [i], e2_);
1225         }
1226
1227         // Closure comparison
1228         BOOST_UBLAS_INLINE
1229         bool same_closure (const vector_binary_scalar2 &vbs2) const {
1230             return (*this).e1_.same_closure (vbs2.e1_) &&
1231                    &e2_ == &(vbs2.e2_);
1232         }
1233
1234         // Iterator types
1235     private:
1236         typedef typename expression1_type::const_iterator const_subiterator1_type;
1237         typedef expression2_type const_subiterator2_type;
1238         typedef const value_type *const_pointer;
1239
1240     public:
1241 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1242         typedef indexed_const_iterator<const_closure_type, typename const_subiterator2_type::iterator_category> const_iterator;
1243         typedef const_iterator iterator;
1244 #else
1245         class const_iterator;
1246         typedef const_iterator iterator;
1247 #endif
1248
1249         // Element lookup
1250         BOOST_UBLAS_INLINE
1251         const_iterator find (size_type i) const {
1252 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1253             const_subiterator1_type it (e1_.find (i));
1254             return const_iterator (*this, it.index ());
1255 #else
1256             return const_iterator (*this, e1_.find (i), const_subiterator2_type (e2_));
1257 #endif
1258         }
1259
1260         // Iterator enhances the iterator of the referenced vector expression
1261         // with the binary functor.
1262
1263 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1264         class const_iterator:
1265             public container_const_reference<vector_binary_scalar2>,
1266             public iterator_base_traits<typename E1::const_iterator::iterator_category>::template
1267                         iterator_base<const_iterator, value_type>::type {
1268         public:
1269             typedef typename E1::const_iterator::iterator_category iterator_category;
1270             typedef typename vector_binary_scalar2::difference_type difference_type;
1271             typedef typename vector_binary_scalar2::value_type value_type;
1272             typedef typename vector_binary_scalar2::const_reference reference;
1273             typedef typename vector_binary_scalar2::const_pointer pointer;
1274
1275             // Construction and destruction
1276             BOOST_UBLAS_INLINE
1277             const_iterator ():
1278                 container_const_reference<self_type> (), it1_ (), it2_ () {}
1279             BOOST_UBLAS_INLINE
1280             const_iterator (const self_type &vbs, const const_subiterator1_type &it1, const const_subiterator2_type &it2):
1281                 container_const_reference<self_type> (vbs), it1_ (it1), it2_ (it2) {}
1282
1283             // Arithmetic
1284             BOOST_UBLAS_INLINE
1285             const_iterator &operator ++ () {
1286                 ++ it1_;
1287                 return *this;
1288             }
1289             BOOST_UBLAS_INLINE
1290             const_iterator &operator -- () {
1291                 -- it1_;
1292                 return *this;
1293             }
1294             BOOST_UBLAS_INLINE
1295             const_iterator &operator += (difference_type n) {
1296                 it1_ += n;
1297                 return *this;
1298             }
1299             BOOST_UBLAS_INLINE
1300             const_iterator &operator -= (difference_type n) {
1301                 it1_ -= n;
1302                 return *this;
1303             }
1304             BOOST_UBLAS_INLINE
1305             difference_type operator - (const const_iterator &it) const {
1306                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1307                 // FIXME we shouldn't compare floats
1308                 // BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
1309                 return it1_ - it.it1_;
1310             }
1311
1312             // Dereference
1313             BOOST_UBLAS_INLINE
1314             const_reference operator * () const {
1315                 return functor_type::apply (*it1_, it2_);
1316             }
1317             BOOST_UBLAS_INLINE
1318             const_reference operator [] (difference_type n) const {
1319                 return *(*this + n);
1320             }
1321
1322             // Index
1323             BOOST_UBLAS_INLINE
1324             size_type index () const {
1325                 return it1_.index ();
1326             }
1327
1328             // Assignment
1329             BOOST_UBLAS_INLINE
1330             const_iterator &operator = (const const_iterator &it) {
1331                 container_const_reference<self_type>::assign (&it ());
1332                 it1_ = it.it1_;
1333                 it2_ = it.it2_;
1334                 return *this;
1335             }
1336
1337             // Comparison
1338             BOOST_UBLAS_INLINE
1339             bool operator == (const const_iterator &it) const {
1340                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1341                 // FIXME we shouldn't compare floats
1342                 // BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
1343                 return it1_ == it.it1_;
1344             }
1345             BOOST_UBLAS_INLINE
1346             bool operator < (const const_iterator &it) const {
1347                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1348                 // FIXME we shouldn't compare floats
1349                 // BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
1350                 return it1_ < it.it1_;
1351             }
1352
1353         private:
1354             const_subiterator1_type it1_;
1355             const_subiterator2_type it2_;
1356         };
1357 #endif
1358
1359         BOOST_UBLAS_INLINE
1360         const_iterator begin () const {
1361             return find (0);
1362         }
1363         BOOST_UBLAS_INLINE
1364         const_iterator end () const {
1365             return find (size ());
1366         }
1367
1368         // Reverse iterator
1369         typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
1370
1371         BOOST_UBLAS_INLINE
1372         const_reverse_iterator rbegin () const {
1373             return const_reverse_iterator (end ());
1374         }
1375         BOOST_UBLAS_INLINE
1376         const_reverse_iterator rend () const {
1377             return const_reverse_iterator (begin ());
1378         }
1379
1380     private:
1381         expression1_closure_type e1_;
1382         expression2_closure_type e2_;
1383     };
1384
1385     template<class E1, class E2, class F>
1386     struct vector_binary_scalar2_traits {
1387         typedef vector_binary_scalar2<E1, E2, F> expression_type;   // allow E2 to be builtin type
1388 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
1389         typedef expression_type result_type;
1390 #else
1391         typedef typename E1::vector_temporary_type result_type;
1392 #endif
1393     };
1394
1395     // (v * t) [i] = v [i] * t
1396     template<class E1, class T2>
1397     BOOST_UBLAS_INLINE
1398     typename vector_binary_scalar2_traits<E1, const T2, scalar_multiplies<typename E1::value_type, T2> >::result_type
1399     operator * (const vector_expression<E1> &e1,
1400                 const T2 &e2) {
1401         typedef typename vector_binary_scalar2_traits<E1, const T2, scalar_multiplies<typename E1::value_type, T2> >::expression_type expression_type;
1402         return expression_type (e1 (), e2);
1403     }
1404
1405     // (v / t) [i] = v [i] / t
1406     template<class E1, class T2>
1407     BOOST_UBLAS_INLINE
1408     typename vector_binary_scalar2_traits<E1, const T2, scalar_divides<typename E1::value_type, T2> >::result_type
1409     operator / (const vector_expression<E1> &e1,
1410                 const T2 &e2) {
1411         typedef typename vector_binary_scalar2_traits<E1, const T2, scalar_divides<typename E1::value_type, T2> >::expression_type expression_type;
1412         return expression_type (e1 (), e2);
1413     }
1414
1415
1416     template<class E, class F>
1417     class vector_scalar_unary:
1418         public scalar_expression<vector_scalar_unary<E, F> > {
1419
1420         typedef E expression_type;
1421         typedef F functor_type;
1422         typedef typename E::const_closure_type expression_closure_type;
1423         typedef typename E::const_iterator::iterator_category iterator_category;
1424         typedef vector_scalar_unary<E, F> self_type;
1425     public:
1426         typedef typename F::result_type value_type;
1427         typedef typename E::difference_type difference_type;
1428         typedef const self_type const_closure_type;
1429         typedef const_closure_type closure_type;
1430         typedef unknown_storage_tag storage_category;
1431
1432         // Construction and destruction
1433         BOOST_UBLAS_INLINE
1434         explicit vector_scalar_unary (const expression_type &e):
1435             e_ (e) {}
1436
1437     private:
1438         // Expression accessors
1439         BOOST_UBLAS_INLINE
1440         const expression_closure_type &expression () const {
1441             return e_;
1442         }
1443
1444     public:
1445         BOOST_UBLAS_INLINE
1446         operator value_type () const {
1447             return evaluate (iterator_category ());
1448         }
1449
1450     private:
1451         // Dense random access specialization
1452         BOOST_UBLAS_INLINE
1453         value_type evaluate (dense_random_access_iterator_tag) const {
1454 #ifdef BOOST_UBLAS_USE_INDEXING
1455             return functor_type::apply (e_);
1456 #elif BOOST_UBLAS_USE_ITERATING
1457             difference_type size = e_.size ();
1458             return functor_type::apply (size, e_.begin ());
1459 #else
1460             difference_type size = e_.size ();
1461             if (size >= BOOST_UBLAS_ITERATOR_THRESHOLD)
1462                 return functor_type::apply (size, e_.begin ());
1463             else
1464                 return functor_type::apply (e_);
1465 #endif
1466         }
1467
1468         // Packed bidirectional specialization
1469         BOOST_UBLAS_INLINE
1470         value_type evaluate (packed_random_access_iterator_tag) const {
1471             return functor_type::apply (e_.begin (), e_.end ());
1472         }
1473
1474         // Sparse bidirectional specialization
1475         BOOST_UBLAS_INLINE
1476         value_type evaluate (sparse_bidirectional_iterator_tag) const {
1477             return functor_type::apply (e_.begin (), e_.end ());
1478         }
1479
1480     private:
1481         expression_closure_type e_;
1482     };
1483
1484     template<class E, class F>
1485     struct vector_scalar_unary_traits {
1486         typedef vector_scalar_unary<E, F> expression_type;
1487 #if !defined (BOOST_UBLAS_SIMPLE_ET_DEBUG) && defined (BOOST_UBLAS_USE_SCALAR_ET)
1488 // FIXME don't define USE_SCALAR_ET other then for testing
1489 // They do not work for complex types
1490          typedef expression_type result_type;
1491 #else
1492          typedef typename F::result_type result_type;
1493 #endif
1494     };
1495
1496     // sum v = sum (v [i])
1497     template<class E>
1498     BOOST_UBLAS_INLINE
1499     typename vector_scalar_unary_traits<E, vector_sum<E> >::result_type
1500     sum (const vector_expression<E> &e) {
1501         typedef typename vector_scalar_unary_traits<E, vector_sum<E> >::expression_type expression_type;
1502         return expression_type (e ());
1503     }
1504
1505     // real: norm_1 v = sum (abs (v [i]))
1506     // complex: norm_1 v = sum (abs (real (v [i])) + abs (imag (v [i])))
1507     template<class E>
1508     BOOST_UBLAS_INLINE
1509     typename vector_scalar_unary_traits<E, vector_norm_1<E> >::result_type
1510     norm_1 (const vector_expression<E> &e) {
1511         typedef typename vector_scalar_unary_traits<E, vector_norm_1<E> >::expression_type expression_type;
1512         return expression_type (e ());
1513     }
1514
1515     // real: norm_2 v = sqrt (sum (v [i] * v [i]))
1516     // complex: norm_2 v = sqrt (sum (v [i] * conj (v [i])))
1517     template<class E>
1518     BOOST_UBLAS_INLINE
1519     typename vector_scalar_unary_traits<E, vector_norm_2<E> >::result_type
1520     norm_2 (const vector_expression<E> &e) {
1521         typedef typename vector_scalar_unary_traits<E, vector_norm_2<E> >::expression_type expression_type;
1522         return expression_type (e ());
1523     }
1524
1525     // real: norm_inf v = maximum (abs (v [i]))
1526     // complex: norm_inf v = maximum (maximum (abs (real (v [i])), abs (imag (v [i]))))
1527     template<class E>
1528     BOOST_UBLAS_INLINE
1529     typename vector_scalar_unary_traits<E, vector_norm_inf<E> >::result_type
1530     norm_inf (const vector_expression<E> &e) {
1531         typedef typename vector_scalar_unary_traits<E, vector_norm_inf<E> >::expression_type expression_type;
1532         return expression_type (e ());
1533     }
1534
1535     // real: index_norm_inf v = minimum (i: abs (v [i]) == maximum (abs (v [i])))
1536     template<class E>
1537     BOOST_UBLAS_INLINE
1538     typename vector_scalar_unary_traits<E, vector_index_norm_inf<E> >::result_type
1539     index_norm_inf (const vector_expression<E> &e) {
1540         typedef typename vector_scalar_unary_traits<E, vector_index_norm_inf<E> >::expression_type expression_type;
1541         return expression_type (e ());
1542     }
1543
1544     template<class E1, class E2, class F>
1545     class vector_scalar_binary:
1546         public scalar_expression<vector_scalar_binary<E1, E2, F> > {
1547
1548         typedef E1 expression1_type;
1549         typedef E2 expression2_type;
1550         typedef F functor_type;
1551         typedef typename E1::const_closure_type expression1_closure_type;
1552         typedef typename E2::const_closure_type expression2_closure_type;
1553         typedef typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
1554                                                   typename E2::const_iterator::iterator_category>::iterator_category iterator_category;
1555         typedef vector_scalar_binary<E1, E2, F> self_type;
1556     public:
1557         static const unsigned complexity = 1;
1558         typedef typename F::result_type value_type;
1559         typedef typename E1::difference_type difference_type;
1560         typedef const self_type const_closure_type;
1561         typedef const_closure_type closure_type;
1562         typedef unknown_storage_tag storage_category;
1563
1564         // Construction and destruction
1565         BOOST_UBLAS_INLINE
1566         vector_scalar_binary (const expression1_type &e1, const expression2_type  &e2):
1567             e1_ (e1), e2_ (e2) {}
1568
1569     private:
1570         // Accessors
1571         BOOST_UBLAS_INLINE
1572         const expression1_closure_type &expression1 () const {
1573             return e1_;
1574         }
1575         BOOST_UBLAS_INLINE
1576         const expression2_closure_type &expression2 () const {
1577             return e2_;
1578         }
1579
1580     public:
1581         BOOST_UBLAS_INLINE
1582         operator value_type () const {
1583             return evaluate (iterator_category ());
1584         }
1585
1586     private:
1587         // Dense random access specialization
1588         BOOST_UBLAS_INLINE
1589         value_type evaluate (dense_random_access_iterator_tag) const {
1590             BOOST_UBLAS_CHECK (e1_.size () == e2_.size (), external_logic());
1591 #ifdef BOOST_UBLAS_USE_INDEXING
1592             return functor_type::apply (e1_, e2_);
1593 #elif BOOST_UBLAS_USE_ITERATING
1594             difference_type size = BOOST_UBLAS_SAME (e1_.size (), e2_.size ());
1595             return functor_type::apply (size, e1_.begin (), e2_.begin ());
1596 #else
1597             difference_type size = BOOST_UBLAS_SAME (e1_.size (), e2_.size ());
1598             if (size >= BOOST_UBLAS_ITERATOR_THRESHOLD)
1599                 return functor_type::apply (size, e1_.begin (), e2_.begin ());
1600             else
1601                 return functor_type::apply (e1_, e2_);
1602 #endif
1603         }
1604
1605         // Packed bidirectional specialization
1606         BOOST_UBLAS_INLINE
1607         value_type evaluate (packed_random_access_iterator_tag) const {
1608             BOOST_UBLAS_CHECK (e1_.size () == e2_.size (), external_logic());
1609             return functor_type::apply (e1_.begin (), e1_.end (), e2_.begin (), e2_.end ());
1610         }
1611
1612         // Sparse bidirectional specialization
1613         BOOST_UBLAS_INLINE
1614         value_type evaluate (sparse_bidirectional_iterator_tag) const {
1615             BOOST_UBLAS_CHECK (e1_.size () == e2_.size (), external_logic());
1616             return functor_type::apply (e1_.begin (), e1_.end (), e2_.begin (), e2_.end (), sparse_bidirectional_iterator_tag ());
1617         }
1618
1619     private:
1620         expression1_closure_type e1_;
1621         expression2_closure_type e2_;
1622     };
1623
1624     template<class E1, class E2, class F>
1625     struct vector_scalar_binary_traits {
1626         typedef vector_scalar_binary<E1, E2, F> expression_type;
1627 #if !defined (BOOST_UBLAS_SIMPLE_ET_DEBUG) && defined (BOOST_UBLAS_USE_SCALAR_ET)
1628 // FIXME don't define USE_SCALAR_ET other then for testing
1629 // They do not work for complex types
1630         typedef expression_type result_type;
1631 #else
1632         typedef typename F::result_type result_type;
1633 #endif
1634     };
1635
1636     // inner_prod (v1, v2) = sum (v1 [i] * v2 [i])
1637     template<class E1, class E2>
1638     BOOST_UBLAS_INLINE
1639     typename vector_scalar_binary_traits<E1, E2, vector_inner_prod<E1, E2,
1640                                                                    typename promote_traits<typename E1::value_type,
1641                                                                                            typename E2::value_type>::promote_type> >::result_type
1642     inner_prod (const vector_expression<E1> &e1,
1643                 const vector_expression<E2> &e2) {
1644         typedef typename vector_scalar_binary_traits<E1, E2, vector_inner_prod<E1, E2,
1645                                                                    typename promote_traits<typename E1::value_type,
1646                                                                                            typename E2::value_type>::promote_type> >::expression_type expression_type;
1647         return expression_type (e1 (), e2 ());
1648     }
1649
1650     template<class E1, class E2>
1651     BOOST_UBLAS_INLINE
1652     typename vector_scalar_binary_traits<E1, E2, vector_inner_prod<E1, E2,
1653                                                                    typename type_traits<typename promote_traits<typename E1::value_type,
1654                                                                                                                 typename E2::value_type>::promote_type>::precision_type> >::result_type
1655     prec_inner_prod (const vector_expression<E1> &e1,
1656                      const vector_expression<E2> &e2) {
1657         typedef typename vector_scalar_binary_traits<E1, E2, vector_inner_prod<E1, E2,
1658                                                                    typename type_traits<typename promote_traits<typename E1::value_type,
1659                                                                                                                 typename E2::value_type>::promote_type>::precision_type> >::expression_type expression_type;
1660         return expression_type (e1 (), e2 ());
1661     }
1662
1663 }}}
1664
1665 #endif