Boost.uBlas 1.49
Linear Algebra in C++: matrices, vectors and numeric algorithms

vector_expression.hpp

Go to the documentation of this file.
00001 //
00002 //  Copyright (c) 2000-2002
00003 //  Joerg Walter, Mathias Koch
00004 //
00005 //  Distributed under the Boost Software License, Version 1.0. (See
00006 //  accompanying file LICENSE_1_0.txt or copy at
00007 //  http://www.boost.org/LICENSE_1_0.txt)
00008 //
00009 //  The authors gratefully acknowledge the support of
00010 //  GeNeSys mbH & Co. KG in producing this work.
00011 //
00012 
00013 #ifndef _BOOST_UBLAS_VECTOR_EXPRESSION_
00014 #define _BOOST_UBLAS_VECTOR_EXPRESSION_
00015 
00016 #include <boost/numeric/ublas/expression_types.hpp>
00017 
00018 
00019 // Expression templates based on ideas of Todd Veldhuizen and Geoffrey Furnish
00020 // Iterators based on ideas of Jeremy Siek
00021 //
00022 // Classes that model the Vector Expression concept
00023 
00024 namespace boost { namespace numeric { namespace ublas {
00025 
00026     template<class E>
00027     class vector_reference:
00028         public vector_expression<vector_reference<E> > {
00029 
00030         typedef vector_reference<E> self_type;
00031     public:
00032 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
00033         using vector_expression<vector_reference<E> >::operator ();
00034 #endif
00035         typedef typename E::size_type size_type;
00036         typedef typename E::difference_type difference_type;
00037         typedef typename E::value_type value_type;
00038         typedef typename E::const_reference const_reference;
00039         typedef typename boost::mpl::if_<boost::is_const<E>,
00040                                           typename E::const_reference,
00041                                           typename E::reference>::type reference;
00042         typedef E referred_type;
00043         typedef const self_type const_closure_type;
00044         typedef self_type closure_type;
00045         typedef typename E::storage_category storage_category;
00046 
00047         // Construction and destruction
00048         BOOST_UBLAS_INLINE
00049         explicit vector_reference (referred_type &e):
00050             e_ (e) {}
00051 
00052         // Accessors
00053         BOOST_UBLAS_INLINE
00054         size_type size () const {
00055             return expression ().size ();
00056         }
00057 
00058     public:
00059         // Expression accessors - const correct
00060         BOOST_UBLAS_INLINE
00061         const referred_type &expression () const {
00062             return e_;
00063         }
00064         BOOST_UBLAS_INLINE
00065         referred_type &expression () {
00066             return e_;
00067         }
00068 
00069     public:
00070         // Element access
00071 #ifndef BOOST_UBLAS_REFERENCE_CONST_MEMBER
00072         BOOST_UBLAS_INLINE
00073         const_reference operator () (size_type i) const {
00074             return expression () (i);
00075         }
00076         BOOST_UBLAS_INLINE
00077         reference operator () (size_type i) {
00078             return expression () (i);
00079         }
00080 
00081         BOOST_UBLAS_INLINE
00082         const_reference operator [] (size_type i) const {
00083             return expression () [i];
00084         }
00085         BOOST_UBLAS_INLINE
00086         reference operator [] (size_type i) {
00087             return expression () [i];
00088         }
00089 #else
00090         BOOST_UBLAS_INLINE
00091         reference operator () (size_type i) const {
00092             return expression () (i);
00093         }
00094 
00095         BOOST_UBLAS_INLINE
00096         reference operator [] (size_type i) const {
00097             return expression () [i];
00098         }
00099 #endif
00100 
00101         // Assignment
00102         BOOST_UBLAS_INLINE
00103         vector_reference &operator = (const vector_reference &v) {
00104             expression ().operator = (v);
00105             return *this;
00106         }
00107         template<class AE>
00108         BOOST_UBLAS_INLINE
00109         vector_reference &operator = (const vector_expression<AE> &ae) {
00110             expression ().operator = (ae);
00111             return *this;
00112         }
00113         template<class AE>
00114         BOOST_UBLAS_INLINE
00115         vector_reference &assign (const vector_expression<AE> &ae) {
00116             expression ().assign (ae);
00117             return *this;
00118         }
00119         template<class AE>
00120         BOOST_UBLAS_INLINE
00121         vector_reference &operator += (const vector_expression<AE> &ae) {
00122             expression ().operator += (ae);
00123             return *this;
00124         }
00125         template<class AE>
00126         BOOST_UBLAS_INLINE
00127         vector_reference &plus_assign (const vector_expression<AE> &ae) {
00128             expression ().plus_assign (ae);
00129             return *this;
00130         }
00131         template<class AE>
00132         BOOST_UBLAS_INLINE
00133         vector_reference &operator -= (const vector_expression<AE> &ae) {
00134             expression ().operator -= (ae);
00135             return *this;
00136         }
00137         template<class AE>
00138         BOOST_UBLAS_INLINE
00139         vector_reference &minus_assign (const vector_expression<AE> &ae) {
00140             expression ().minus_assign (ae);
00141             return *this;
00142         }
00143         template<class AT>
00144         BOOST_UBLAS_INLINE
00145         vector_reference &operator *= (const AT &at) {
00146             expression ().operator *= (at);
00147             return *this;
00148         }
00149         template<class AT>
00150         BOOST_UBLAS_INLINE
00151         vector_reference &operator /= (const AT &at) {
00152             expression ().operator /= (at);
00153             return *this;
00154         }
00155 
00156         // Swapping
00157         BOOST_UBLAS_INLINE
00158         void swap (vector_reference &v) {
00159             expression ().swap (v.expression ());
00160         }
00161 
00162         // Closure comparison
00163         BOOST_UBLAS_INLINE
00164         bool same_closure (const vector_reference &vr) const {
00165             return &(*this).e_ == &vr.e_;
00166         }
00167 
00168         // Iterator types
00169         typedef typename E::const_iterator const_iterator;
00170         typedef typename boost::mpl::if_<boost::is_const<E>,
00171                                           typename E::const_iterator,
00172                                           typename E::iterator>::type iterator;
00173 
00174         // Element lookup
00175         BOOST_UBLAS_INLINE
00176         const_iterator find (size_type i) const {
00177             return expression ().find (i);
00178         }
00179         BOOST_UBLAS_INLINE
00180         iterator find (size_type i) {
00181             return expression ().find (i);
00182         }
00183 
00184         // Iterator is the iterator of the referenced expression.
00185 
00186         BOOST_UBLAS_INLINE
00187         const_iterator begin () const {
00188             return expression ().begin ();
00189         }
00190         BOOST_UBLAS_INLINE
00191         const_iterator end () const {
00192             return expression ().end ();
00193         }
00194 
00195         BOOST_UBLAS_INLINE
00196         iterator begin () {
00197             return expression ().begin ();
00198         }
00199         BOOST_UBLAS_INLINE
00200         iterator end () {
00201             return expression ().end ();
00202         }
00203 
00204         // Reverse iterator
00205         typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
00206         typedef reverse_iterator_base<iterator> reverse_iterator;
00207 
00208         BOOST_UBLAS_INLINE
00209         const_reverse_iterator rbegin () const {
00210             return const_reverse_iterator (end ());
00211         }
00212         BOOST_UBLAS_INLINE
00213         const_reverse_iterator rend () const {
00214             return const_reverse_iterator (begin ());
00215         }
00216         BOOST_UBLAS_INLINE
00217         reverse_iterator rbegin () {
00218             return reverse_iterator (end ());
00219         }
00220         BOOST_UBLAS_INLINE
00221         reverse_iterator rend () {
00222             return reverse_iterator (begin ());
00223         }
00224 
00225     private:
00226         referred_type &e_;
00227     };
00228 
00229 
00230     template<class E, class F>
00231     class vector_unary:
00232         public vector_expression<vector_unary<E, F> > {
00233 
00234         typedef F functor_type;
00235         typedef typename boost::mpl::if_<boost::is_same<F, scalar_identity<typename E::value_type> >,
00236                                           E,
00237                                           const E>::type expression_type;
00238         typedef typename boost::mpl::if_<boost::is_const<expression_type>,
00239                                           typename E::const_closure_type,
00240                                           typename E::closure_type>::type expression_closure_type;
00241         typedef vector_unary<E, F> self_type;
00242     public:
00243 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
00244         using vector_expression<vector_unary<E, F> >::operator ();
00245 #endif
00246         typedef typename E::size_type size_type;
00247         typedef typename E::difference_type difference_type;
00248         typedef typename F::result_type value_type;
00249         typedef value_type const_reference;
00250         typedef typename boost::mpl::if_<boost::is_same<F, scalar_identity<value_type> >,
00251                                           typename E::reference,
00252                                           value_type>::type reference;
00253         typedef const self_type const_closure_type;
00254         typedef self_type closure_type;
00255         typedef unknown_storage_tag storage_category;
00256 
00257         // Construction and destruction
00258         BOOST_UBLAS_INLINE
00259         // May be used as mutable expression.
00260         explicit vector_unary (expression_type &e):
00261             e_ (e) {}
00262 
00263         // Accessors
00264         BOOST_UBLAS_INLINE
00265         size_type size () const {
00266             return e_.size ();
00267         }
00268 
00269     public:
00270         // Expression accessors
00271         BOOST_UBLAS_INLINE
00272         const expression_closure_type &expression () const {
00273             return e_;
00274         }
00275 
00276     public:
00277         // Element access
00278         BOOST_UBLAS_INLINE
00279         const_reference operator () (size_type i) const {
00280             return functor_type::apply (e_ (i));
00281         }
00282         BOOST_UBLAS_INLINE
00283         reference operator () (size_type i) {
00284             BOOST_STATIC_ASSERT ((boost::is_same<functor_type, scalar_identity<value_type > >::value));
00285             return e_ (i);
00286         }
00287 
00288         BOOST_UBLAS_INLINE
00289         const_reference operator [] (size_type i) const {
00290             return functor_type::apply (e_ [i]);
00291         }
00292         BOOST_UBLAS_INLINE
00293         reference operator [] (size_type i) {
00294             BOOST_STATIC_ASSERT ((boost::is_same<functor_type, scalar_identity<value_type > >::value));
00295             return e_ [i];
00296         }
00297 
00298         // Closure comparison
00299         BOOST_UBLAS_INLINE
00300         bool same_closure (const vector_unary &vu) const {
00301             return (*this).expression ().same_closure (vu.expression ());
00302         }
00303 
00304         // Iterator types
00305     private:
00306         typedef typename E::const_iterator const_subiterator_type;
00307         typedef const value_type *const_pointer;
00308 
00309     public:
00310 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
00311         typedef indexed_const_iterator<const_closure_type, typename const_subiterator_type::iterator_category> const_iterator;
00312         typedef const_iterator iterator;
00313 #else
00314         class const_iterator;
00315         typedef const_iterator iterator;
00316 #endif
00317 
00318         // Element lookup
00319         BOOST_UBLAS_INLINE
00320         const_iterator find (size_type i) const {
00321 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
00322             const_subiterator_type it (e_.find (i));
00323             return const_iterator (*this, it.index ());
00324 #else
00325             return const_iterator (*this, e_.find (i));
00326 #endif
00327         }
00328 
00329         // Iterator enhances the iterator of the referenced expression
00330         // with the unary functor.
00331 
00332 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
00333         class const_iterator:
00334             public container_const_reference<vector_unary>,
00335             public iterator_base_traits<typename E::const_iterator::iterator_category>::template
00336                         iterator_base<const_iterator, value_type>::type {
00337         public:
00338             typedef typename E::const_iterator::iterator_category iterator_category;
00339             typedef typename vector_unary::difference_type difference_type;
00340             typedef typename vector_unary::value_type value_type;
00341             typedef typename vector_unary::const_reference reference;
00342             typedef typename vector_unary::const_pointer pointer;
00343 
00344             // Construction and destruction
00345             BOOST_UBLAS_INLINE
00346             const_iterator ():
00347                 container_const_reference<self_type> (), it_ () {}
00348             BOOST_UBLAS_INLINE
00349             const_iterator (const self_type &vu, const const_subiterator_type &it):
00350                 container_const_reference<self_type> (vu), it_ (it) {}
00351 
00352             // Arithmetic
00353             BOOST_UBLAS_INLINE
00354             const_iterator &operator ++ () {
00355                 ++ it_;
00356                 return *this;
00357             }
00358             BOOST_UBLAS_INLINE
00359             const_iterator &operator -- () {
00360                 -- it_;
00361                 return *this;
00362             }
00363             BOOST_UBLAS_INLINE
00364             const_iterator &operator += (difference_type n) {
00365                 it_ += n;
00366                 return *this;
00367             }
00368             BOOST_UBLAS_INLINE
00369             const_iterator &operator -= (difference_type n) {
00370                 it_ -= n;
00371                 return *this;
00372             }
00373             BOOST_UBLAS_INLINE
00374             difference_type operator - (const const_iterator &it) const {
00375                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
00376                 return it_ - it.it_;
00377             }
00378 
00379             // Dereference
00380             BOOST_UBLAS_INLINE
00381             const_reference operator * () const {
00382                 return functor_type::apply (*it_);
00383             }
00384             BOOST_UBLAS_INLINE
00385             const_reference operator [] (difference_type n) const {
00386                 return *(*this + n);
00387             }
00388 
00389             // Index
00390             BOOST_UBLAS_INLINE
00391             size_type index () const {
00392                 return it_.index ();
00393             }
00394 
00395             // Assignment
00396             BOOST_UBLAS_INLINE
00397             const_iterator &operator = (const const_iterator &it) {
00398                 container_const_reference<self_type>::assign (&it ());
00399                 it_ = it.it_;
00400                 return *this;
00401             }
00402 
00403             // Comparison
00404             BOOST_UBLAS_INLINE
00405             bool operator == (const const_iterator &it) const {
00406                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
00407                 return it_ == it.it_;
00408             }
00409             BOOST_UBLAS_INLINE
00410             bool operator < (const const_iterator &it) const {
00411                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
00412                 return it_ < it.it_;
00413             }
00414 
00415         private:
00416             const_subiterator_type it_;
00417         };
00418 #endif
00419 
00420         BOOST_UBLAS_INLINE
00421         const_iterator begin () const {
00422             return find (0); 
00423         }
00424         BOOST_UBLAS_INLINE
00425         const_iterator end () const {
00426             return find (size ());
00427         }
00428 
00429         // Reverse iterator
00430         typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
00431 
00432         BOOST_UBLAS_INLINE
00433         const_reverse_iterator rbegin () const {
00434             return const_reverse_iterator (end ());
00435         }
00436         BOOST_UBLAS_INLINE
00437         const_reverse_iterator rend () const {
00438             return const_reverse_iterator (begin ());
00439         }
00440 
00441     private:
00442         expression_closure_type e_;
00443     };
00444 
00445     template<class E, class F>
00446     struct vector_unary_traits {
00447         typedef vector_unary<E, F> expression_type;
00448 //FIXME
00449 // #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
00450         typedef expression_type result_type;
00451 // #else
00452 //         typedef typename E::vector_temporary_type result_type;
00453 // #endif
00454     };
00455 
00456     // (- v) [i] = - v [i]
00457     template<class E> 
00458     BOOST_UBLAS_INLINE
00459     typename vector_unary_traits<E, scalar_negate<typename E::value_type> >::result_type
00460     operator - (const vector_expression<E> &e) {
00461         typedef typename vector_unary_traits<E, scalar_negate<typename E::value_type> >::expression_type expression_type;
00462         return expression_type (e ());
00463     }
00464 
00465     // (conj v) [i] = conj (v [i])
00466     template<class E> 
00467     BOOST_UBLAS_INLINE
00468     typename vector_unary_traits<E, scalar_conj<typename E::value_type> >::result_type
00469     conj (const vector_expression<E> &e) {
00470         typedef typename vector_unary_traits<E, scalar_conj<typename E::value_type> >::expression_type expression_type;
00471         return expression_type (e ());
00472     }
00473 
00474     // (real v) [i] = real (v [i])
00475     template<class E>
00476     BOOST_UBLAS_INLINE
00477     typename vector_unary_traits<E, scalar_real<typename E::value_type> >::result_type
00478     real (const vector_expression<E> &e) {
00479         typedef typename vector_unary_traits<E, scalar_real<typename E::value_type> >::expression_type expression_type;
00480         return expression_type (e ());
00481     }
00482 
00483     // (imag v) [i] = imag (v [i])
00484     template<class E>
00485     BOOST_UBLAS_INLINE
00486     typename vector_unary_traits<E, scalar_imag<typename E::value_type> >::result_type
00487     imag (const vector_expression<E> &e) {
00488         typedef typename vector_unary_traits<E, scalar_imag<typename E::value_type> >::expression_type expression_type;
00489         return expression_type (e ());
00490     }
00491 
00492     // (trans v) [i] = v [i]
00493     template<class E>
00494     BOOST_UBLAS_INLINE
00495     typename vector_unary_traits<const E, scalar_identity<typename E::value_type> >::result_type
00496     trans (const vector_expression<E> &e) {
00497         typedef typename vector_unary_traits<const E, scalar_identity<typename E::value_type> >::expression_type expression_type;
00498         return expression_type (e ());
00499     }
00500     template<class E>
00501     BOOST_UBLAS_INLINE
00502     typename vector_unary_traits<E, scalar_identity<typename E::value_type> >::result_type
00503     trans (vector_expression<E> &e) {
00504         typedef typename vector_unary_traits<E, scalar_identity<typename E::value_type> >::expression_type expression_type;
00505         return expression_type (e ());
00506     }
00507 
00508     // (herm v) [i] = conj (v [i])
00509     template<class E>
00510     BOOST_UBLAS_INLINE
00511     typename vector_unary_traits<E, scalar_conj<typename E::value_type> >::result_type
00512     herm (const vector_expression<E> &e) {
00513         typedef typename vector_unary_traits<E, scalar_conj<typename E::value_type> >::expression_type expression_type;
00514         return expression_type (e ());
00515     }
00516 
00517     template<class E1, class E2, class F>
00518     class vector_binary:
00519         public vector_expression<vector_binary<E1, E2, F> > {
00520 
00521         typedef E1 expression1_type;
00522         typedef E2 expression2_type;
00523         typedef F functor_type;
00524         typedef typename E1::const_closure_type expression1_closure_type;
00525         typedef typename E2::const_closure_type expression2_closure_type;
00526         typedef vector_binary<E1, E2, F> self_type;
00527     public:
00528 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
00529         using vector_expression<vector_binary<E1, E2, F> >::operator ();
00530 #endif
00531         typedef typename promote_traits<typename E1::size_type, typename E2::size_type>::promote_type size_type;
00532         typedef typename promote_traits<typename E1::difference_type, typename E2::difference_type>::promote_type difference_type;
00533         typedef typename F::result_type value_type;
00534         typedef value_type const_reference;
00535         typedef const_reference reference;
00536         typedef const self_type const_closure_type;
00537         typedef const_closure_type closure_type;
00538         typedef unknown_storage_tag storage_category;
00539 
00540         // Construction and destruction
00541         BOOST_UBLAS_INLINE
00542         vector_binary (const expression1_type &e1, const expression2_type &e2):
00543             e1_ (e1), e2_ (e2) {}
00544 
00545         // Accessors
00546         BOOST_UBLAS_INLINE
00547         size_type size () const { 
00548             return BOOST_UBLAS_SAME (e1_.size (), e2_.size ()); 
00549         }
00550 
00551     private:
00552         // Accessors
00553         BOOST_UBLAS_INLINE
00554         const expression1_closure_type &expression1 () const {
00555             return e1_;
00556         }
00557         BOOST_UBLAS_INLINE
00558         const expression2_closure_type &expression2 () const {
00559             return e2_;
00560         }
00561 
00562     public:
00563         // Element access
00564         BOOST_UBLAS_INLINE
00565         const_reference operator () (size_type i) const {
00566             return functor_type::apply (e1_ (i), e2_ (i));
00567         }
00568 
00569         BOOST_UBLAS_INLINE
00570         const_reference operator [] (size_type i) const {
00571             return functor_type::apply (e1_ [i], e2_ [i]);
00572         }
00573 
00574         // Closure comparison
00575         BOOST_UBLAS_INLINE
00576         bool same_closure (const vector_binary &vb) const {
00577             return (*this).expression1 ().same_closure (vb.expression1 ()) &&
00578                    (*this).expression2 ().same_closure (vb.expression2 ());
00579         }
00580 
00581         // Iterator types
00582     private:
00583         typedef typename E1::const_iterator const_subiterator1_type;
00584         typedef typename E2::const_iterator const_subiterator2_type;
00585         typedef const value_type *const_pointer;
00586 
00587     public:
00588 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
00589         typedef typename iterator_restrict_traits<typename const_subiterator1_type::iterator_category,
00590                                                   typename const_subiterator2_type::iterator_category>::iterator_category iterator_category;
00591         typedef indexed_const_iterator<const_closure_type, iterator_category> const_iterator;
00592         typedef const_iterator iterator;
00593 #else
00594         class const_iterator;
00595         typedef const_iterator iterator;
00596 #endif
00597 
00598         // Element lookup
00599         BOOST_UBLAS_INLINE
00600         const_iterator find (size_type i) const {
00601             const_subiterator1_type it1 (e1_.find (i));
00602             const_subiterator1_type it1_end (e1_.find (size ()));
00603             const_subiterator2_type it2 (e2_.find (i));
00604             const_subiterator2_type it2_end (e2_.find (size ()));
00605             i = (std::min) (it1 != it1_end ? it1.index () : size (),
00606                           it2 != it2_end ? it2.index () : size ());
00607 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
00608             return const_iterator (*this, i);
00609 #else
00610             return const_iterator (*this, i, it1, it1_end, it2, it2_end);
00611 #endif
00612         }
00613 
00614         // Iterator merges the iterators of the referenced expressions and
00615         // enhances them with the binary functor.
00616 
00617 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
00618         class const_iterator:
00619             public container_const_reference<vector_binary>,
00620             public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
00621                                                                           typename E2::const_iterator::iterator_category>::iterator_category>::template
00622                         iterator_base<const_iterator, value_type>::type {
00623         public:
00624             typedef typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
00625                                                       typename E2::const_iterator::iterator_category>::iterator_category iterator_category;
00626             typedef typename vector_binary::difference_type difference_type;
00627             typedef typename vector_binary::value_type value_type;
00628             typedef typename vector_binary::const_reference reference;
00629             typedef typename vector_binary::const_pointer pointer;
00630 
00631             // Construction and destruction
00632             BOOST_UBLAS_INLINE
00633             const_iterator ():
00634                 container_const_reference<self_type> (), i_ (), it1_ (), it1_end_ (), it2_ (), it2_end_ () {}
00635             BOOST_UBLAS_INLINE
00636             const_iterator (const self_type &vb, size_type i,
00637                             const const_subiterator1_type &it1, const const_subiterator1_type &it1_end,
00638                             const const_subiterator2_type &it2, const const_subiterator2_type &it2_end):
00639                 container_const_reference<self_type> (vb), i_ (i), it1_ (it1), it1_end_ (it1_end), it2_ (it2), it2_end_ (it2_end) {}
00640 
00641         private: 
00642             // Dense specializations
00643             BOOST_UBLAS_INLINE
00644             void increment (dense_random_access_iterator_tag) {
00645                 ++ i_; ++ it1_; ++ it2_;
00646             }
00647             BOOST_UBLAS_INLINE
00648             void decrement (dense_random_access_iterator_tag) {
00649                 -- i_; -- it1_; -- it2_;
00650             }
00651             BOOST_UBLAS_INLINE
00652             void increment (dense_random_access_iterator_tag, difference_type n) {
00653                 i_ += n; it1_ += n; it2_ += n;
00654             }
00655             BOOST_UBLAS_INLINE
00656             void decrement (dense_random_access_iterator_tag, difference_type n) {
00657                 i_ -= n; it1_ -= n; it2_ -= n;
00658             }
00659             BOOST_UBLAS_INLINE
00660             value_type dereference (dense_random_access_iterator_tag) const {
00661                 return functor_type::apply (*it1_, *it2_);
00662             }
00663 
00664             // Packed specializations
00665             BOOST_UBLAS_INLINE
00666             void increment (packed_random_access_iterator_tag) {
00667                 if (it1_ != it1_end_)
00668                     if (it1_.index () <= i_)
00669                         ++ it1_;
00670                 if (it2_ != it2_end_)
00671                     if (it2_.index () <= i_)
00672                         ++ it2_;
00673                 ++ i_;
00674             }
00675             BOOST_UBLAS_INLINE
00676             void decrement (packed_random_access_iterator_tag) {
00677                 if (it1_ != it1_end_)
00678                     if (i_ <= it1_.index ())
00679                         -- it1_;
00680                 if (it2_ != it2_end_)
00681                     if (i_ <= it2_.index ())
00682                         -- it2_;
00683                 -- i_;
00684             }
00685             BOOST_UBLAS_INLINE
00686             void increment (packed_random_access_iterator_tag, difference_type n) {
00687                 while (n > 0) {
00688                     increment (packed_random_access_iterator_tag ());
00689                     --n;
00690                 }
00691                 while (n < 0) {
00692                     decrement (packed_random_access_iterator_tag ());
00693                     ++n;
00694                 }
00695             }
00696             BOOST_UBLAS_INLINE
00697             void decrement (packed_random_access_iterator_tag, difference_type n) {
00698                 while (n > 0) {
00699                     decrement (packed_random_access_iterator_tag ());
00700                     --n;
00701                 }
00702                 while (n < 0) {
00703                     increment (packed_random_access_iterator_tag ());
00704                     ++n;
00705                 }
00706             }
00707             BOOST_UBLAS_INLINE
00708             value_type dereference (packed_random_access_iterator_tag) const {
00709                 value_type t1 = value_type/*zero*/();
00710                 if (it1_ != it1_end_)
00711                     if (it1_.index () == i_)
00712                         t1 = *it1_;
00713                 value_type t2 = value_type/*zero*/();
00714                 if (it2_ != it2_end_)
00715                     if (it2_.index () == i_)
00716                         t2 = *it2_;
00717                 return functor_type::apply (t1, t2);
00718             }
00719 
00720             // Sparse specializations
00721             BOOST_UBLAS_INLINE
00722             void increment (sparse_bidirectional_iterator_tag) {
00723                 size_type index1 = (*this) ().size ();
00724                 if (it1_ != it1_end_) {
00725                     if  (it1_.index () <= i_)
00726                         ++ it1_;
00727                     if (it1_ != it1_end_)
00728                         index1 = it1_.index ();
00729                 }
00730                 size_type index2 = (*this) ().size ();
00731                 if (it2_ != it2_end_) {
00732                     if (it2_.index () <= i_)
00733                         ++ it2_;
00734                     if (it2_ != it2_end_)
00735                         index2 = it2_.index ();
00736                 }
00737                 i_ = (std::min) (index1, index2);
00738             }
00739             BOOST_UBLAS_INLINE
00740             void decrement (sparse_bidirectional_iterator_tag) {
00741                 size_type index1 = (*this) ().size ();
00742                 if (it1_ != it1_end_) {
00743                     if (i_ <= it1_.index ())
00744                         -- it1_;
00745                     if (it1_ != it1_end_)
00746                         index1 = it1_.index ();
00747                 }
00748                 size_type index2 = (*this) ().size ();
00749                 if (it2_ != it2_end_) {
00750                     if (i_ <= it2_.index ())
00751                         -- it2_;
00752                     if (it2_ != it2_end_)
00753                         index2 = it2_.index ();
00754                 }
00755                 i_ = (std::max) (index1, index2);
00756             }
00757             BOOST_UBLAS_INLINE
00758             void increment (sparse_bidirectional_iterator_tag, difference_type n) {
00759                 while (n > 0) {
00760                     increment (sparse_bidirectional_iterator_tag ());
00761                     --n;
00762                 }
00763                 while (n < 0) {
00764                     decrement (sparse_bidirectional_iterator_tag ());
00765                     ++n;
00766                 }
00767             }
00768             BOOST_UBLAS_INLINE
00769             void decrement (sparse_bidirectional_iterator_tag, difference_type n) {
00770                 while (n > 0) {
00771                     decrement (sparse_bidirectional_iterator_tag ());
00772                     --n;
00773                 }
00774                 while (n < 0) {
00775                     increment (sparse_bidirectional_iterator_tag ());
00776                     ++n;
00777                 }
00778             }
00779             BOOST_UBLAS_INLINE
00780             value_type dereference (sparse_bidirectional_iterator_tag) const {
00781                 value_type t1 = value_type/*zero*/();
00782                 if (it1_ != it1_end_)
00783                     if (it1_.index () == i_)
00784                         t1 = *it1_;
00785                 value_type t2 = value_type/*zero*/();
00786                 if (it2_ != it2_end_)
00787                     if (it2_.index () == i_)
00788                         t2 = *it2_;
00789                 return functor_type::apply (t1, t2);
00790             }
00791 
00792         public: 
00793             // Arithmetic
00794             BOOST_UBLAS_INLINE
00795             const_iterator &operator ++ () {
00796                 increment (iterator_category ());
00797                 return *this;
00798             }
00799             BOOST_UBLAS_INLINE
00800             const_iterator &operator -- () {
00801                 decrement (iterator_category ());
00802                 return *this;
00803             }
00804             BOOST_UBLAS_INLINE
00805             const_iterator &operator += (difference_type n) {
00806                 increment (iterator_category (), n);
00807                 return *this;
00808             }
00809             BOOST_UBLAS_INLINE
00810             const_iterator &operator -= (difference_type n) {
00811                 decrement (iterator_category (), n);
00812                 return *this;
00813             }
00814             BOOST_UBLAS_INLINE
00815             difference_type operator - (const const_iterator &it) const {
00816                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
00817                 return index () - it.index ();
00818             }
00819 
00820             // Dereference
00821             BOOST_UBLAS_INLINE
00822             const_reference operator * () const {
00823                 return dereference (iterator_category ());
00824             }
00825             BOOST_UBLAS_INLINE
00826             const_reference operator [] (difference_type n) const {
00827                 return *(*this + n);
00828             }
00829 
00830             // Index
00831             BOOST_UBLAS_INLINE
00832             size_type index () const {
00833                 return i_;
00834             }
00835 
00836             // Assignment
00837             BOOST_UBLAS_INLINE
00838             const_iterator &operator = (const const_iterator &it) {
00839                 container_const_reference<self_type>::assign (&it ());
00840                 i_ = it.i_;
00841                 it1_ = it.it1_;
00842                 it1_end_ = it.it1_end_;
00843                 it2_ = it.it2_;
00844                 it2_end_ = it.it2_end_;
00845                 return *this;
00846             }
00847 
00848             // Comparison
00849             BOOST_UBLAS_INLINE
00850             bool operator == (const const_iterator &it) const {
00851                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
00852                 return index () == it.index ();
00853             }
00854             BOOST_UBLAS_INLINE
00855             bool operator < (const const_iterator &it) const {
00856                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
00857                 return index () < it.index ();
00858             }
00859 
00860         private:
00861             size_type i_;
00862             const_subiterator1_type it1_;
00863             const_subiterator1_type it1_end_;
00864             const_subiterator2_type it2_;
00865             const_subiterator2_type it2_end_;
00866         };
00867 #endif
00868 
00869         BOOST_UBLAS_INLINE
00870         const_iterator begin () const {
00871             return find (0);
00872         }
00873         BOOST_UBLAS_INLINE
00874         const_iterator end () const {
00875             return find (size ());
00876         }
00877 
00878         // Reverse iterator
00879         typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
00880 
00881         BOOST_UBLAS_INLINE
00882         const_reverse_iterator rbegin () const {
00883             return const_reverse_iterator (end ());
00884         }
00885         BOOST_UBLAS_INLINE
00886         const_reverse_iterator rend () const {
00887             return const_reverse_iterator (begin ());
00888         }
00889 
00890     private:
00891         expression1_closure_type e1_;
00892         expression2_closure_type e2_;
00893     };
00894 
00895     template<class E1, class E2, class F>
00896     struct vector_binary_traits {
00897         typedef vector_binary<E1, E2, F> expression_type;
00898 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
00899         typedef expression_type result_type; 
00900 #else
00901         typedef typename E1::vector_temporary_type result_type;
00902 #endif
00903     };
00904 
00905     // (v1 + v2) [i] = v1 [i] + v2 [i]
00906     template<class E1, class E2>
00907     BOOST_UBLAS_INLINE
00908     typename vector_binary_traits<E1, E2, scalar_plus<typename E1::value_type, 
00909                                                       typename E2::value_type> >::result_type
00910     operator + (const vector_expression<E1> &e1,
00911                 const vector_expression<E2> &e2) {
00912         typedef typename vector_binary_traits<E1, E2, scalar_plus<typename E1::value_type,
00913                                                                               typename E2::value_type> >::expression_type expression_type;
00914         return expression_type (e1 (), e2 ());
00915     }
00916 
00917     // (v1 - v2) [i] = v1 [i] - v2 [i]
00918     template<class E1, class E2>
00919     BOOST_UBLAS_INLINE
00920     typename vector_binary_traits<E1, E2, scalar_minus<typename E1::value_type,
00921                                                        typename E2::value_type> >::result_type
00922     operator - (const vector_expression<E1> &e1,
00923                 const vector_expression<E2> &e2) {
00924         typedef typename vector_binary_traits<E1, E2, scalar_minus<typename E1::value_type,
00925                                                                                typename E2::value_type> >::expression_type expression_type;
00926         return expression_type (e1 (), e2 ());
00927     }
00928 
00929     // (v1 * v2) [i] = v1 [i] * v2 [i]
00930     template<class E1, class E2>
00931     BOOST_UBLAS_INLINE
00932     typename vector_binary_traits<E1, E2, scalar_multiplies<typename E1::value_type,
00933                                                             typename E2::value_type> >::result_type
00934     element_prod (const vector_expression<E1> &e1,
00935                   const vector_expression<E2> &e2) {
00936         typedef typename vector_binary_traits<E1, E2, scalar_multiplies<typename E1::value_type,
00937                                                                                     typename E2::value_type> >::expression_type expression_type;
00938         return expression_type (e1 (), e2 ());
00939     }
00940 
00941     // (v1 / v2) [i] = v1 [i] / v2 [i]
00942     template<class E1, class E2>
00943     BOOST_UBLAS_INLINE
00944     typename vector_binary_traits<E1, E2, scalar_divides<typename E1::value_type,
00945                                                          typename E2::value_type> >::result_type
00946     element_div (const vector_expression<E1> &e1,
00947                  const vector_expression<E2> &e2) {
00948         typedef typename vector_binary_traits<E1, E2, scalar_divides<typename E1::value_type,
00949                                                                                  typename E2::value_type> >::expression_type expression_type;
00950         return expression_type (e1 (), e2 ());
00951     }
00952 
00953 
00954     template<class E1, class E2, class F>
00955     class vector_binary_scalar1:
00956         public vector_expression<vector_binary_scalar1<E1, E2, F> > {
00957 
00958         typedef F functor_type;
00959         typedef E1 expression1_type;
00960         typedef E2 expression2_type;
00961     public:
00962         typedef const E1& expression1_closure_type;
00963         typedef typename E2::const_closure_type expression2_closure_type;
00964     private:
00965         typedef vector_binary_scalar1<E1, E2, F> self_type;
00966     public:
00967 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
00968         using vector_expression<vector_binary_scalar1<E1, E2, F> >::operator ();
00969 #endif
00970         typedef typename E2::size_type size_type;
00971         typedef typename E2::difference_type difference_type;
00972         typedef typename F::result_type value_type;
00973         typedef value_type const_reference;
00974         typedef const_reference reference;
00975         typedef const self_type const_closure_type;
00976         typedef const_closure_type closure_type;
00977         typedef unknown_storage_tag storage_category;
00978 
00979         // Construction and destruction
00980         BOOST_UBLAS_INLINE
00981         vector_binary_scalar1 (const expression1_type &e1, const expression2_type &e2):
00982             e1_ (e1), e2_ (e2) {}
00983 
00984         // Accessors
00985         BOOST_UBLAS_INLINE
00986         size_type size () const {
00987             return e2_.size ();
00988         }
00989 
00990     public:
00991         // Element access
00992         BOOST_UBLAS_INLINE
00993         const_reference operator () (size_type i) const {
00994             return functor_type::apply (e1_, e2_ (i));
00995         }
00996 
00997         BOOST_UBLAS_INLINE
00998         const_reference operator [] (size_type i) const {
00999             return functor_type::apply (e1_, e2_ [i]);
01000         }
01001 
01002         // Closure comparison
01003         BOOST_UBLAS_INLINE
01004         bool same_closure (const vector_binary_scalar1 &vbs1) const {
01005             return &e1_ == &(vbs1.e1_) &&
01006                    (*this).e2_.same_closure (vbs1.e2_);
01007         }
01008 
01009         // Iterator types
01010     private:
01011         typedef expression1_type const_subiterator1_type;
01012         typedef typename expression2_type::const_iterator const_subiterator2_type;
01013         typedef const value_type *const_pointer;
01014 
01015     public:
01016 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
01017         typedef indexed_const_iterator<const_closure_type, typename const_subiterator2_type::iterator_category> const_iterator;
01018         typedef const_iterator iterator;
01019 #else
01020         class const_iterator;
01021         typedef const_iterator iterator;
01022 #endif
01023 
01024         // Element lookup
01025         BOOST_UBLAS_INLINE
01026         const_iterator find (size_type i) const {
01027 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
01028             const_subiterator2_type it (e2_.find (i));
01029             return const_iterator (*this, it.index ());
01030 #else
01031             return const_iterator (*this, const_subiterator1_type (e1_), e2_.find (i));
01032 #endif
01033         }
01034 
01035         // Iterator enhances the iterator of the referenced vector expression
01036         // with the binary functor.
01037 
01038 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
01039         class const_iterator:
01040             public container_const_reference<vector_binary_scalar1>,
01041             public iterator_base_traits<typename E2::const_iterator::iterator_category>::template
01042                         iterator_base<const_iterator, value_type>::type {
01043         public:
01044             typedef typename E2::const_iterator::iterator_category iterator_category;
01045             typedef typename vector_binary_scalar1::difference_type difference_type;
01046             typedef typename vector_binary_scalar1::value_type value_type;
01047             typedef typename vector_binary_scalar1::const_reference reference;
01048             typedef typename vector_binary_scalar1::const_pointer pointer;
01049 
01050             // Construction and destruction
01051             BOOST_UBLAS_INLINE
01052             const_iterator ():
01053                 container_const_reference<self_type> (), it1_ (), it2_ () {}
01054             BOOST_UBLAS_INLINE
01055             const_iterator (const self_type &vbs, const const_subiterator1_type &it1, const const_subiterator2_type &it2):
01056                 container_const_reference<self_type> (vbs), it1_ (it1), it2_ (it2) {}
01057 
01058             // Arithmetic
01059             BOOST_UBLAS_INLINE
01060             const_iterator &operator ++ () {
01061                 ++ it2_;
01062                 return *this;
01063             }
01064             BOOST_UBLAS_INLINE
01065             const_iterator &operator -- () {
01066                 -- it2_;
01067                 return *this;
01068             }
01069             BOOST_UBLAS_INLINE
01070             const_iterator &operator += (difference_type n) {
01071                 it2_ += n;
01072                 return *this;
01073             }
01074             BOOST_UBLAS_INLINE
01075             const_iterator &operator -= (difference_type n) {
01076                 it2_ -= n;
01077                 return *this;
01078             }
01079             BOOST_UBLAS_INLINE
01080             difference_type operator - (const const_iterator &it) const {
01081                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
01082                 // FIXME we shouldn't compare floats
01083                 // BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
01084                 return it2_ - it.it2_;
01085             }
01086 
01087             // Dereference
01088             BOOST_UBLAS_INLINE
01089             const_reference operator * () const {
01090                 return functor_type::apply (it1_, *it2_);
01091             }
01092             BOOST_UBLAS_INLINE
01093             const_reference operator [] (difference_type n) const {
01094                 return *(*this + n);
01095             }
01096 
01097             // Index
01098             BOOST_UBLAS_INLINE
01099             size_type index () const {
01100                 return it2_.index ();
01101             }
01102 
01103             // Assignment 
01104             BOOST_UBLAS_INLINE
01105             const_iterator &operator = (const const_iterator &it) {
01106                 container_const_reference<self_type>::assign (&it ());
01107                 it1_ = it.it1_;
01108                 it2_ = it.it2_;
01109                 return *this;
01110             }
01111 
01112             // Comparison
01113             BOOST_UBLAS_INLINE
01114             bool operator == (const const_iterator &it) const {
01115                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
01116                 // FIXME we shouldn't compare floats
01117                 // BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
01118                 return it2_ == it.it2_;
01119             }
01120             BOOST_UBLAS_INLINE
01121             bool operator < (const const_iterator &it) const {
01122                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
01123                 // FIXME we shouldn't compare floats
01124                 // BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
01125                 return it2_ < it.it2_;
01126             }
01127 
01128         private:
01129             const_subiterator1_type it1_;
01130             const_subiterator2_type it2_;
01131         };
01132 #endif
01133 
01134         BOOST_UBLAS_INLINE
01135         const_iterator begin () const {
01136             return find (0); 
01137         }
01138         BOOST_UBLAS_INLINE
01139         const_iterator end () const {
01140             return find (size ()); 
01141         }
01142 
01143         // Reverse iterator
01144         typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
01145 
01146         BOOST_UBLAS_INLINE
01147         const_reverse_iterator rbegin () const {
01148             return const_reverse_iterator (end ());
01149         }
01150         BOOST_UBLAS_INLINE
01151         const_reverse_iterator rend () const {
01152             return const_reverse_iterator (begin ());
01153         }
01154 
01155     private:
01156         expression1_closure_type e1_;
01157         expression2_closure_type e2_;
01158     };
01159 
01160     template<class E1, class E2, class F>
01161     struct vector_binary_scalar1_traits {
01162         typedef vector_binary_scalar1<E1, E2, F> expression_type;   // allow E1 to be builtin type
01163 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
01164         typedef expression_type result_type;
01165 #else
01166         typedef typename E2::vector_temporary_type result_type;
01167 #endif
01168     };
01169 
01170     // (t * v) [i] = t * v [i]
01171     template<class T1, class E2>
01172     BOOST_UBLAS_INLINE
01173     typename enable_if< is_convertible<T1, typename E2::value_type >,    
01174     typename vector_binary_scalar1_traits<const T1, E2, scalar_multiplies<T1, typename E2::value_type> >::result_type
01175     >::type
01176     operator * (const T1 &e1,
01177                 const vector_expression<E2> &e2) {
01178         typedef typename vector_binary_scalar1_traits<const T1, E2, scalar_multiplies<T1, typename E2::value_type> >::expression_type expression_type;
01179         return expression_type (e1, e2 ());
01180     }
01181 
01182 
01183     template<class E1, class E2, class F>
01184     class vector_binary_scalar2:
01185         public vector_expression<vector_binary_scalar2<E1, E2, F> > {
01186 
01187         typedef F functor_type;
01188         typedef E1 expression1_type;
01189         typedef E2 expression2_type;
01190         typedef typename E1::const_closure_type expression1_closure_type;
01191         typedef const E2& expression2_closure_type;
01192         typedef vector_binary_scalar2<E1, E2, F> self_type;
01193     public:
01194 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
01195         using vector_expression<vector_binary_scalar2<E1, E2, F> >::operator ();
01196 #endif
01197         typedef typename E1::size_type size_type;
01198         typedef typename E1::difference_type difference_type;
01199         typedef typename F::result_type value_type;
01200         typedef value_type const_reference;
01201         typedef const_reference reference;
01202         typedef const self_type const_closure_type;
01203         typedef const_closure_type closure_type;
01204         typedef unknown_storage_tag storage_category;
01205 
01206         // Construction and destruction
01207         BOOST_UBLAS_INLINE
01208         vector_binary_scalar2 (const expression1_type &e1, const expression2_type &e2):
01209             e1_ (e1), e2_ (e2) {}
01210 
01211         // Accessors
01212         BOOST_UBLAS_INLINE
01213         size_type size () const {
01214             return e1_.size (); 
01215         }
01216 
01217     public:
01218         // Element access
01219         BOOST_UBLAS_INLINE
01220         const_reference operator () (size_type i) const {
01221             return functor_type::apply (e1_ (i), e2_);
01222         }
01223 
01224         BOOST_UBLAS_INLINE
01225         const_reference operator [] (size_type i) const {
01226             return functor_type::apply (e1_ [i], e2_);
01227         }
01228 
01229         // Closure comparison
01230         BOOST_UBLAS_INLINE
01231         bool same_closure (const vector_binary_scalar2 &vbs2) const {
01232             return (*this).e1_.same_closure (vbs2.e1_) &&
01233                    &e2_ == &(vbs2.e2_);
01234         }
01235 
01236         // Iterator types
01237     private:
01238         typedef typename expression1_type::const_iterator const_subiterator1_type;
01239         typedef expression2_type const_subiterator2_type;
01240         typedef const value_type *const_pointer;
01241 
01242     public:
01243 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
01244         typedef indexed_const_iterator<const_closure_type, typename const_subiterator2_type::iterator_category> const_iterator;
01245         typedef const_iterator iterator;
01246 #else
01247         class const_iterator;
01248         typedef const_iterator iterator;
01249 #endif
01250 
01251         // Element lookup
01252         BOOST_UBLAS_INLINE
01253         const_iterator find (size_type i) const {
01254 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
01255             const_subiterator1_type it (e1_.find (i));
01256             return const_iterator (*this, it.index ());
01257 #else
01258             return const_iterator (*this, e1_.find (i), const_subiterator2_type (e2_));
01259 #endif
01260         }
01261 
01262         // Iterator enhances the iterator of the referenced vector expression
01263         // with the binary functor.
01264 
01265 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
01266         class const_iterator:
01267             public container_const_reference<vector_binary_scalar2>,
01268             public iterator_base_traits<typename E1::const_iterator::iterator_category>::template
01269                         iterator_base<const_iterator, value_type>::type {
01270         public:
01271             typedef typename E1::const_iterator::iterator_category iterator_category;
01272             typedef typename vector_binary_scalar2::difference_type difference_type;
01273             typedef typename vector_binary_scalar2::value_type value_type;
01274             typedef typename vector_binary_scalar2::const_reference reference;
01275             typedef typename vector_binary_scalar2::const_pointer pointer;
01276 
01277             // Construction and destruction
01278             BOOST_UBLAS_INLINE
01279             const_iterator ():
01280                 container_const_reference<self_type> (), it1_ (), it2_ () {}
01281             BOOST_UBLAS_INLINE
01282             const_iterator (const self_type &vbs, const const_subiterator1_type &it1, const const_subiterator2_type &it2):
01283                 container_const_reference<self_type> (vbs), it1_ (it1), it2_ (it2) {}
01284 
01285             // Arithmetic
01286             BOOST_UBLAS_INLINE
01287             const_iterator &operator ++ () {
01288                 ++ it1_;
01289                 return *this;
01290             }
01291             BOOST_UBLAS_INLINE
01292             const_iterator &operator -- () {
01293                 -- it1_;
01294                 return *this;
01295             }
01296             BOOST_UBLAS_INLINE
01297             const_iterator &operator += (difference_type n) {
01298                 it1_ += n;
01299                 return *this;
01300             }
01301             BOOST_UBLAS_INLINE
01302             const_iterator &operator -= (difference_type n) {
01303                 it1_ -= n;
01304                 return *this;
01305             }
01306             BOOST_UBLAS_INLINE
01307             difference_type operator - (const const_iterator &it) const {
01308                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
01309                 // FIXME we shouldn't compare floats
01310                 // BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
01311                 return it1_ - it.it1_;
01312             }
01313 
01314             // Dereference
01315             BOOST_UBLAS_INLINE
01316             const_reference operator * () const {
01317                 return functor_type::apply (*it1_, it2_);
01318             }
01319             BOOST_UBLAS_INLINE
01320             const_reference operator [] (difference_type n) const {
01321                 return *(*this + n);
01322             }
01323 
01324             // Index
01325             BOOST_UBLAS_INLINE
01326             size_type index () const {
01327                 return it1_.index ();
01328             }
01329 
01330             // Assignment
01331             BOOST_UBLAS_INLINE
01332             const_iterator &operator = (const const_iterator &it) {
01333                 container_const_reference<self_type>::assign (&it ());
01334                 it1_ = it.it1_;
01335                 it2_ = it.it2_;
01336                 return *this;
01337             }
01338 
01339             // Comparison
01340             BOOST_UBLAS_INLINE
01341             bool operator == (const const_iterator &it) const {
01342                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
01343                 // FIXME we shouldn't compare floats
01344                 // BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
01345                 return it1_ == it.it1_;
01346             }
01347             BOOST_UBLAS_INLINE
01348             bool operator < (const const_iterator &it) const {
01349                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
01350                 // FIXME we shouldn't compare floats
01351                 // BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
01352                 return it1_ < it.it1_;
01353             }
01354 
01355         private:
01356             const_subiterator1_type it1_;
01357             const_subiterator2_type it2_;
01358         };
01359 #endif
01360 
01361         BOOST_UBLAS_INLINE
01362         const_iterator begin () const {
01363             return find (0);
01364         }
01365         BOOST_UBLAS_INLINE
01366         const_iterator end () const {
01367             return find (size ());
01368         }
01369 
01370         // Reverse iterator
01371         typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
01372 
01373         BOOST_UBLAS_INLINE
01374         const_reverse_iterator rbegin () const {
01375             return const_reverse_iterator (end ());
01376         }
01377         BOOST_UBLAS_INLINE
01378         const_reverse_iterator rend () const {
01379             return const_reverse_iterator (begin ());
01380         }
01381 
01382     private:
01383         expression1_closure_type e1_;
01384         expression2_closure_type e2_;
01385     };
01386 
01387     template<class E1, class E2, class F>
01388     struct vector_binary_scalar2_traits {
01389         typedef vector_binary_scalar2<E1, E2, F> expression_type;   // allow E2 to be builtin type
01390 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
01391         typedef expression_type result_type;
01392 #else
01393         typedef typename E1::vector_temporary_type result_type;
01394 #endif
01395     };
01396 
01397     // (v * t) [i] = v [i] * t
01398     template<class E1, class T2>
01399     BOOST_UBLAS_INLINE
01400     typename enable_if< is_convertible<T2, typename E1::value_type >,    
01401     typename vector_binary_scalar2_traits<E1, const T2, scalar_multiplies<typename E1::value_type, T2> >::result_type
01402     >::type
01403     operator * (const vector_expression<E1> &e1,
01404                 const T2 &e2) {
01405         typedef typename vector_binary_scalar2_traits<E1, const T2, scalar_multiplies<typename E1::value_type, T2> >::expression_type expression_type;
01406         return expression_type (e1 (), e2);
01407     }
01408 
01409     // (v / t) [i] = v [i] / t
01410     template<class E1, class T2>
01411     BOOST_UBLAS_INLINE
01412     typename vector_binary_scalar2_traits<E1, const T2, scalar_divides<typename E1::value_type, T2> >::result_type
01413     operator / (const vector_expression<E1> &e1,
01414                 const T2 &e2) {
01415         typedef typename vector_binary_scalar2_traits<E1, const T2, scalar_divides<typename E1::value_type, T2> >::expression_type expression_type;
01416         return expression_type (e1 (), e2);
01417     }
01418 
01419 
01420     template<class E, class F>
01421     class vector_scalar_unary:
01422         public scalar_expression<vector_scalar_unary<E, F> > {
01423 
01424         typedef E expression_type;
01425         typedef F functor_type;
01426         typedef typename E::const_closure_type expression_closure_type;
01427         typedef typename E::const_iterator::iterator_category iterator_category;
01428         typedef vector_scalar_unary<E, F> self_type;
01429     public:
01430         typedef typename F::result_type value_type;
01431         typedef typename E::difference_type difference_type;
01432         typedef const self_type const_closure_type;
01433         typedef const_closure_type closure_type;
01434         typedef unknown_storage_tag storage_category;
01435 
01436         // Construction and destruction
01437         BOOST_UBLAS_INLINE
01438         explicit vector_scalar_unary (const expression_type &e):
01439             e_ (e) {}
01440 
01441     private:
01442         // Expression accessors
01443         BOOST_UBLAS_INLINE
01444         const expression_closure_type &expression () const {
01445             return e_;
01446         }
01447 
01448     public:
01449         BOOST_UBLAS_INLINE
01450         operator value_type () const {
01451             return evaluate (iterator_category ());
01452         }
01453 
01454     private:
01455         // Dense random access specialization
01456         BOOST_UBLAS_INLINE
01457         value_type evaluate (dense_random_access_iterator_tag) const {
01458 #ifdef BOOST_UBLAS_USE_INDEXING
01459             return functor_type::apply (e_);
01460 #elif BOOST_UBLAS_USE_ITERATING
01461             difference_type size = e_.size ();
01462             return functor_type::apply (size, e_.begin ());
01463 #else
01464             difference_type size = e_.size ();
01465             if (size >= BOOST_UBLAS_ITERATOR_THRESHOLD)
01466                 return functor_type::apply (size, e_.begin ());
01467             else
01468                 return functor_type::apply (e_);
01469 #endif
01470         }
01471 
01472         // Packed bidirectional specialization
01473         BOOST_UBLAS_INLINE
01474         value_type evaluate (packed_random_access_iterator_tag) const {
01475             return functor_type::apply (e_.begin (), e_.end ());
01476         }
01477 
01478         // Sparse bidirectional specialization
01479         BOOST_UBLAS_INLINE
01480         value_type evaluate (sparse_bidirectional_iterator_tag) const {
01481             return functor_type::apply (e_.begin (), e_.end ());
01482         }
01483 
01484     private:
01485         expression_closure_type e_;
01486     };
01487 
01488     template<class E, class F>
01489     struct vector_scalar_unary_traits {
01490         typedef vector_scalar_unary<E, F> expression_type;
01491 #if !defined (BOOST_UBLAS_SIMPLE_ET_DEBUG) && defined (BOOST_UBLAS_USE_SCALAR_ET)
01492 // FIXME don't define USE_SCALAR_ET other then for testing
01493 // They do not work for complex types
01494          typedef expression_type result_type;
01495 #else
01496          typedef typename F::result_type result_type;
01497 #endif
01498     };
01499 
01500     // sum v = sum (v [i])
01501     template<class E>
01502     BOOST_UBLAS_INLINE
01503     typename vector_scalar_unary_traits<E, vector_sum<E> >::result_type
01504     sum (const vector_expression<E> &e) {
01505         typedef typename vector_scalar_unary_traits<E, vector_sum<E> >::expression_type expression_type;
01506         return expression_type (e ());
01507     }
01508 
01509     // real: norm_1 v = sum (abs (v [i]))
01510     // complex: norm_1 v = sum (abs (real (v [i])) + abs (imag (v [i])))
01511     template<class E>
01512     BOOST_UBLAS_INLINE
01513     typename vector_scalar_unary_traits<E, vector_norm_1<E> >::result_type
01514     norm_1 (const vector_expression<E> &e) {
01515         typedef typename vector_scalar_unary_traits<E, vector_norm_1<E> >::expression_type expression_type;
01516         return expression_type (e ());
01517     }
01518 
01519     // real: norm_2 v = sqrt (sum (v [i] * v [i]))
01520     // complex: norm_2 v = sqrt (sum (v [i] * conj (v [i])))
01521     template<class E>
01522     BOOST_UBLAS_INLINE
01523     typename vector_scalar_unary_traits<E, vector_norm_2<E> >::result_type
01524     norm_2 (const vector_expression<E> &e) {
01525         typedef typename vector_scalar_unary_traits<E, vector_norm_2<E> >::expression_type expression_type;
01526         return expression_type (e ());
01527     }
01528 
01529     // real: norm_inf v = maximum (abs (v [i]))
01530     // complex: norm_inf v = maximum (maximum (abs (real (v [i])), abs (imag (v [i]))))
01531     template<class E>
01532     BOOST_UBLAS_INLINE
01533     typename vector_scalar_unary_traits<E, vector_norm_inf<E> >::result_type
01534     norm_inf (const vector_expression<E> &e) {
01535         typedef typename vector_scalar_unary_traits<E, vector_norm_inf<E> >::expression_type expression_type;
01536         return expression_type (e ());
01537     }
01538 
01539     // real: index_norm_inf v = minimum (i: abs (v [i]) == maximum (abs (v [i])))
01540     template<class E>
01541     BOOST_UBLAS_INLINE
01542     typename vector_scalar_unary_traits<E, vector_index_norm_inf<E> >::result_type
01543     index_norm_inf (const vector_expression<E> &e) {
01544         typedef typename vector_scalar_unary_traits<E, vector_index_norm_inf<E> >::expression_type expression_type;
01545         return expression_type (e ());
01546     }
01547 
01548     template<class E1, class E2, class F>
01549     class vector_scalar_binary:
01550         public scalar_expression<vector_scalar_binary<E1, E2, F> > {
01551 
01552         typedef E1 expression1_type;
01553         typedef E2 expression2_type;
01554         typedef F functor_type;
01555         typedef typename E1::const_closure_type expression1_closure_type;
01556         typedef typename E2::const_closure_type expression2_closure_type;
01557         typedef typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
01558                                                   typename E2::const_iterator::iterator_category>::iterator_category iterator_category;
01559         typedef vector_scalar_binary<E1, E2, F> self_type;
01560     public:
01561         static const unsigned complexity = 1;
01562         typedef typename F::result_type value_type;
01563         typedef typename E1::difference_type difference_type;
01564         typedef const self_type const_closure_type;
01565         typedef const_closure_type closure_type;
01566         typedef unknown_storage_tag storage_category;
01567 
01568         // Construction and destruction
01569         BOOST_UBLAS_INLINE
01570         vector_scalar_binary (const expression1_type &e1, const expression2_type  &e2):
01571             e1_ (e1), e2_ (e2) {}
01572 
01573     private:
01574         // Accessors
01575         BOOST_UBLAS_INLINE
01576         const expression1_closure_type &expression1 () const {
01577             return e1_;
01578         }
01579         BOOST_UBLAS_INLINE
01580         const expression2_closure_type &expression2 () const {
01581             return e2_;
01582         }
01583 
01584     public:
01585         BOOST_UBLAS_INLINE
01586         operator value_type () const {
01587             return evaluate (iterator_category ());
01588         }
01589 
01590     private:
01591         // Dense random access specialization
01592         BOOST_UBLAS_INLINE
01593         value_type evaluate (dense_random_access_iterator_tag) const {
01594             BOOST_UBLAS_CHECK (e1_.size () == e2_.size (), external_logic());
01595 #ifdef BOOST_UBLAS_USE_INDEXING
01596             return functor_type::apply (e1_, e2_);
01597 #elif BOOST_UBLAS_USE_ITERATING
01598             difference_type size = BOOST_UBLAS_SAME (e1_.size (), e2_.size ());
01599             return functor_type::apply (size, e1_.begin (), e2_.begin ());
01600 #else
01601             difference_type size = BOOST_UBLAS_SAME (e1_.size (), e2_.size ());
01602             if (size >= BOOST_UBLAS_ITERATOR_THRESHOLD)
01603                 return functor_type::apply (size, e1_.begin (), e2_.begin ());
01604             else
01605                 return functor_type::apply (e1_, e2_);
01606 #endif
01607         }
01608 
01609         // Packed bidirectional specialization
01610         BOOST_UBLAS_INLINE
01611         value_type evaluate (packed_random_access_iterator_tag) const {
01612             BOOST_UBLAS_CHECK (e1_.size () == e2_.size (), external_logic());
01613             return functor_type::apply (e1_.begin (), e1_.end (), e2_.begin (), e2_.end ());
01614         }
01615 
01616         // Sparse bidirectional specialization
01617         BOOST_UBLAS_INLINE
01618         value_type evaluate (sparse_bidirectional_iterator_tag) const {
01619             BOOST_UBLAS_CHECK (e1_.size () == e2_.size (), external_logic());
01620             return functor_type::apply (e1_.begin (), e1_.end (), e2_.begin (), e2_.end (), sparse_bidirectional_iterator_tag ());
01621         }
01622 
01623     private:
01624         expression1_closure_type e1_;
01625         expression2_closure_type e2_;
01626     };
01627 
01628     template<class E1, class E2, class F>
01629     struct vector_scalar_binary_traits {
01630         typedef vector_scalar_binary<E1, E2, F> expression_type;
01631 #if !defined (BOOST_UBLAS_SIMPLE_ET_DEBUG) && defined (BOOST_UBLAS_USE_SCALAR_ET)
01632 // FIXME don't define USE_SCALAR_ET other then for testing
01633 // They do not work for complex types
01634         typedef expression_type result_type;
01635 #else
01636         typedef typename F::result_type result_type;
01637 #endif
01638     };
01639 
01640     // inner_prod (v1, v2) = sum (v1 [i] * v2 [i])
01641     template<class E1, class E2>
01642     BOOST_UBLAS_INLINE
01643     typename vector_scalar_binary_traits<E1, E2, vector_inner_prod<E1, E2,
01644                                                                    typename promote_traits<typename E1::value_type,
01645                                                                                            typename E2::value_type>::promote_type> >::result_type
01646     inner_prod (const vector_expression<E1> &e1,
01647                 const vector_expression<E2> &e2) {
01648         typedef typename vector_scalar_binary_traits<E1, E2, vector_inner_prod<E1, E2,
01649                                                                    typename promote_traits<typename E1::value_type,
01650                                                                                            typename E2::value_type>::promote_type> >::expression_type expression_type;
01651         return expression_type (e1 (), e2 ());
01652     }
01653 
01654     template<class E1, class E2>
01655     BOOST_UBLAS_INLINE
01656     typename vector_scalar_binary_traits<E1, E2, vector_inner_prod<E1, E2,
01657                                                                    typename type_traits<typename promote_traits<typename E1::value_type,
01658                                                                                                                 typename E2::value_type>::promote_type>::precision_type> >::result_type
01659     prec_inner_prod (const vector_expression<E1> &e1,
01660                      const vector_expression<E2> &e2) {
01661         typedef typename vector_scalar_binary_traits<E1, E2, vector_inner_prod<E1, E2,
01662                                                                    typename type_traits<typename promote_traits<typename E1::value_type,
01663                                                                                                                 typename E2::value_type>::promote_type>::precision_type> >::expression_type expression_type;
01664         return expression_type (e1 (), e2 ());
01665     }
01666 
01667 }}}
01668 
01669 #endif