numeric/ublas/matrix_expression.hpp

0001 //
0002 //  Copyright (c) 2000-2002
0003 //  Joerg Walter, Mathias Koch
0004 //
0005 //  Distributed under the Boost Software License, Version 1.0. (See
0006 //  accompanying file LICENSE_1_0.txt or copy at
0007 //  http://www.boost.org/LICENSE_1_0.txt)
0008 //
0009 //  The authors gratefully acknowledge the support of
0010 //  GeNeSys mbH & Co. KG in producing this work.
0011 //
0012
0013 #ifndef _BOOST_UBLAS_MATRIX_EXPRESSION_
0014 #define _BOOST_UBLAS_MATRIX_EXPRESSION_
0015
0016 #include <boost/numeric/ublas/vector_expression.hpp>
0017
0018 // Expression templates based on ideas of Todd Veldhuizen and Geoffrey Furnish
0019 // Iterators based on ideas of Jeremy Siek
0020 //
0021 // Classes that model the Matrix Expression concept
0022
0023 namespace boost { namespace numeric { namespace ublas {
0024
0025 template<class E>
0026 class matrix_reference:
0027     public matrix_expression<matrix_reference<E> > {
0028
0029     typedef matrix_reference<E> self_type;
0030 public:
0031 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
0032     using matrix_expression<self_type>::operator ();
0033 #endif
0034     typedef typename E::size_type size_type;
0035     typedef typename E::difference_type difference_type;
0036     typedef typename E::value_type value_type;
0037     typedef typename E::const_reference const_reference;
0038     typedef typename boost::mpl::if_<boost::is_const<E>,
0039     typename E::const_reference,
0040     typename E::reference>::type reference;
0041     typedef E referred_type;
0042     typedef const self_type const_closure_type;
0043     typedef self_type closure_type;
0044     typedef typename E::orientation_category orientation_category;
0045     typedef typename E::storage_category storage_category;
0046
0047     // Construction and destruction
0048     BOOST_UBLAS_INLINE
0049     explicit matrix_reference (referred_type &e):
0050       e_ (e) {}
0051
0052     // Accessors
0053     BOOST_UBLAS_INLINE
0054     size_type size1 () const {
0055         return e_.size1 ();
0056     }
0057     BOOST_UBLAS_INLINE
0058     size_type size2 () const {
0059         return e_.size2 ();
0060     }
0061
0062 public:
0063     // Expression accessors - const correct
0064     BOOST_UBLAS_INLINE
0065     const referred_type &expression () const {
0066         return e_;
0067     }
0068     BOOST_UBLAS_INLINE
0069     referred_type &expression () {
0070         return e_;
0071     }
0072
0073 public:
0074     // Element access
0075 #ifndef BOOST_UBLAS_REFERENCE_CONST_MEMBER
0076     BOOST_UBLAS_INLINE
0077     const_reference operator () (size_type i, size_type j) const {
0078         return expression () (i, j);
0079     }
0080     BOOST_UBLAS_INLINE
0081     reference operator () (size_type i, size_type j) {
0082         return expression () (i, j);
0083     }
0084 #else
0085     BOOST_UBLAS_INLINE
0086     reference operator () (size_type i, size_type j) const {
0087         return expression () (i, j);
0088     }
0089 #endif
0090
0091 #ifndef BOOST_UBLAS_REFERENCE_CONST_MEMBER
0092     BOOST_UBLAS_INLINE
0093     const_reference operator () (size_type i) const {
0094         return expression () (i);
0095     }
0096     BOOST_UBLAS_INLINE
0097     reference operator () (size_type i) {
0098         return expression () (i);
0099     }
0100 #else
0101     BOOST_UBLAS_INLINE
0102     reference operator () (size_type i) const {
0103         return expression () (i);
0104     }
0105 #endif
0106
0107
0108     // Assignment
0109     BOOST_UBLAS_INLINE
0110     matrix_reference &operator = (const matrix_reference &m) {
0111         expression ().operator = (m);
0112         return *this;
0113     }
0114     template<class AE>
0115     BOOST_UBLAS_INLINE
0116     matrix_reference &operator = (const matrix_expression<AE> &ae) {
0117         expression ().operator = (ae);
0118         return *this;
0119     }
0120     template<class AE>
0121     BOOST_UBLAS_INLINE
0122     matrix_reference &assign (const matrix_expression<AE> &ae) {
0123         expression ().assign (ae);
0124         return *this;
0125     }
0126     template<class AE>
0127     BOOST_UBLAS_INLINE
0128     matrix_reference &operator += (const matrix_expression<AE> &ae) {
0129         expression ().operator += (ae);
0130         return *this;
0131     }
0132     template<class AE>
0133     BOOST_UBLAS_INLINE
0134     matrix_reference &plus_assign (const matrix_expression<AE> &ae) {
0135         expression ().plus_assign (ae);
0136         return *this;
0137     }
0138     template<class AE>
0139     BOOST_UBLAS_INLINE
0140     matrix_reference &operator -= (const matrix_expression<AE> &ae) {
0141         expression ().operator -= (ae);
0142         return *this;
0143     }
0144     template<class AE>
0145     BOOST_UBLAS_INLINE
0146     matrix_reference &minus_assign (const matrix_expression<AE> &ae) {
0147         expression ().minus_assign (ae);
0148         return *this;
0149     }
0150     template<class AT>
0151     BOOST_UBLAS_INLINE
0152     matrix_reference &operator *= (const AT &at) {
0153         expression ().operator *= (at);
0154         return *this;
0155     }
0156     template<class AT>
0157     BOOST_UBLAS_INLINE
0158     matrix_reference &operator /= (const AT &at) {
0159         expression ().operator /= (at);
0160         return *this;
0161     }
0162
0163     // Swapping
0164     BOOST_UBLAS_INLINE
0165     void swap (matrix_reference &m) {
0166         expression ().swap (m.expression ());
0167     }
0168
0169     // Closure comparison
0170     BOOST_UBLAS_INLINE
0171     bool same_closure (const matrix_reference &mr) const {
0172         return &(*this).e_ == &mr.e_;
0173     }
0174
0175     // Iterator types
0176     typedef typename E::const_iterator1 const_iterator1;
0177     typedef typename boost::mpl::if_<boost::is_const<E>,
0178     typename E::const_iterator1,
0179     typename E::iterator1>::type iterator1;
0180     typedef typename E::const_iterator2 const_iterator2;
0181     typedef typename boost::mpl::if_<boost::is_const<E>,
0182     typename E::const_iterator2,
0183     typename E::iterator2>::type iterator2;
0184
0185     // Element lookup
0186     BOOST_UBLAS_INLINE
0187     const_iterator1 find1 (int rank, size_type i, size_type j) const {
0188         return expression ().find1 (rank, i, j);
0189     }
0190     BOOST_UBLAS_INLINE
0191     iterator1 find1 (int rank, size_type i, size_type j) {
0192         return expression ().find1 (rank, i, j);
0193     }
0194     BOOST_UBLAS_INLINE
0195     const_iterator2 find2 (int rank, size_type i, size_type j) const {
0196         return expression ().find2 (rank, i, j);
0197     }
0198     BOOST_UBLAS_INLINE
0199     iterator2 find2 (int rank, size_type i, size_type j) {
0200         return expression ().find2 (rank, i, j);
0201     }
0202
0203     // Iterators are the iterators of the referenced expression.
0204
0205     BOOST_UBLAS_INLINE
0206     const_iterator1 begin1 () const {
0207         return expression ().begin1 ();
0208     }
0209     BOOST_UBLAS_INLINE
0210     const_iterator1 cbegin1 () const {
0211         return begin1 ();
0212     }
0213     BOOST_UBLAS_INLINE
0214     const_iterator1 end1 () const {
0215         return expression ().end1 ();
0216     }
0217     BOOST_UBLAS_INLINE
0218     const_iterator1 cend1 () const {
0219         return end1 ();
0220     }
0221
0222     BOOST_UBLAS_INLINE
0223     iterator1 begin1 () {
0224         return expression ().begin1 ();
0225     }
0226     BOOST_UBLAS_INLINE
0227     iterator1 end1 () {
0228         return expression ().end1 ();
0229     }
0230
0231     BOOST_UBLAS_INLINE
0232     const_iterator2 begin2 () const {
0233         return expression ().begin2 ();
0234     }
0235     BOOST_UBLAS_INLINE
0236     const_iterator2 cbegin2 () const {
0237         return begin2 ();
0238     }
0239     BOOST_UBLAS_INLINE
0240     const_iterator2 end2 () const {
0241         return expression ().end2 ();
0242     }
0243     BOOST_UBLAS_INLINE
0244     const_iterator2 cend2 () const {
0245         return end2 ();
0246     }
0247
0248     BOOST_UBLAS_INLINE
0249     iterator2 begin2 () {
0250         return expression ().begin2 ();
0251     }
0252     BOOST_UBLAS_INLINE
0253     iterator2 end2 () {
0254         return expression ().end2 ();
0255     }
0256
0257     // Reverse iterators
0258     typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
0259     typedef reverse_iterator_base1<iterator1> reverse_iterator1;
0260
0261     BOOST_UBLAS_INLINE
0262     const_reverse_iterator1 rbegin1 () const {
0263         return const_reverse_iterator1 (end1 ());
0264     }
0265     BOOST_UBLAS_INLINE
0266     const_reverse_iterator1 crbegin1 () const {
0267         return rbegin1 ();
0268     }
0269     BOOST_UBLAS_INLINE
0270     const_reverse_iterator1 rend1 () const {
0271         return const_reverse_iterator1 (begin1 ());
0272     }
0273     BOOST_UBLAS_INLINE
0274     const_reverse_iterator1 crend1 () const {
0275         return rend1 ();
0276     }
0277
0278     BOOST_UBLAS_INLINE
0279     reverse_iterator1 rbegin1 () {
0280         return reverse_iterator1 (end1 ());
0281     }
0282     BOOST_UBLAS_INLINE
0283     reverse_iterator1 rend1 () {
0284         return reverse_iterator1 (begin1 ());
0285     }
0286
0287     typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
0288     typedef reverse_iterator_base2<iterator2> reverse_iterator2;
0289
0290     BOOST_UBLAS_INLINE
0291     const_reverse_iterator2 rbegin2 () const {
0292         return const_reverse_iterator2 (end2 ());
0293     }
0294     BOOST_UBLAS_INLINE
0295     const_reverse_iterator2 crbegin2 () const {
0296         return rbegin2 ();
0297     }
0298     BOOST_UBLAS_INLINE
0299     const_reverse_iterator2 rend2 () const {
0300         return const_reverse_iterator2 (begin2 ());
0301     }
0302     BOOST_UBLAS_INLINE
0303     const_reverse_iterator2 crend2 () const {
0304         return rend2 ();
0305     }
0306
0307     BOOST_UBLAS_INLINE
0308     reverse_iterator2 rbegin2 () {
0309         return reverse_iterator2 (end2 ());
0310     }
0311     BOOST_UBLAS_INLINE
0312     reverse_iterator2 rend2 () {
0313         return reverse_iterator2 (begin2 ());
0314     }
0315
0316 private:
0317     referred_type &e_;
0318 };
0319
0320
0321 template<class E1, class E2, class F>
0322 class vector_matrix_binary:
0323     public matrix_expression<vector_matrix_binary<E1, E2, F> > {
0324
0325     typedef E1 expression1_type;
0326     typedef E2 expression2_type;
0327 public:
0328     typedef typename E1::const_closure_type expression1_closure_type;
0329     typedef typename E2::const_closure_type expression2_closure_type;
0330 private:
0331     typedef vector_matrix_binary<E1, E2, F> self_type;
0332 public:
0333 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
0334     using matrix_expression<self_type>::operator ();
0335 #endif
0336     typedef F functor_type;
0337     typedef typename promote_traits<typename E1::size_type, typename E2::size_type>::promote_type size_type;
0338     typedef typename promote_traits<typename E1::difference_type, typename E2::difference_type>::promote_type difference_type;
0339     typedef typename F::result_type value_type;
0340     typedef value_type const_reference;
0341     typedef const_reference reference;
0342     typedef const self_type const_closure_type;
0343     typedef const_closure_type closure_type;
0344     typedef unknown_orientation_tag orientation_category;
0345     typedef unknown_storage_tag storage_category;
0346
0347     // Construction and destruction
0348     BOOST_UBLAS_INLINE
0349     vector_matrix_binary (const expression1_type &e1, const expression2_type &e2):
0350       e1_ (e1), e2_ (e2) {}
0351
0352     // Accessors
0353     BOOST_UBLAS_INLINE
0354     size_type size1 () const {
0355         return e1_.size ();
0356     }
0357     BOOST_UBLAS_INLINE
0358     size_type size2 () const {
0359         return e2_.size ();
0360     }
0361
0362 public:
0363     // Expression accessors
0364     BOOST_UBLAS_INLINE
0365     const expression1_closure_type &expression1 () const {
0366         return e1_;
0367     }
0368     BOOST_UBLAS_INLINE
0369     const expression2_closure_type &expression2 () const {
0370         return e2_;
0371     }
0372
0373 public:
0374     // Element access
0375     BOOST_UBLAS_INLINE
0376     const_reference operator () (size_type i, size_type j) const {
0377         return functor_type::apply (e1_ (i), e2_ (j));
0378     }
0379
0380
0381
0382     // Closure comparison
0383     BOOST_UBLAS_INLINE
0384     bool same_closure (const vector_matrix_binary &vmb) const {
0385         return (*this).expression1 ().same_closure (vmb.expression1 ()) &&
0386             (*this).expression2 ().same_closure (vmb.expression2 ());
0387     }
0388
0389     // Iterator types
0390 private:
0391     typedef typename E1::const_iterator const_subiterator1_type;
0392     typedef typename E2::const_iterator const_subiterator2_type;
0393     typedef const value_type *const_pointer;
0394
0395 public:
0396 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
0397     typedef typename iterator_restrict_traits<typename const_subiterator1_type::iterator_category,
0398     typename const_subiterator2_type::iterator_category>::iterator_category iterator_category;
0399     typedef indexed_const_iterator1<const_closure_type, iterator_category> const_iterator1;
0400     typedef const_iterator1 iterator1;
0401     typedef indexed_const_iterator2<const_closure_type, iterator_category> const_iterator2;
0402     typedef const_iterator2 iterator2;
0403 #else
0404     class const_iterator1;
0405     typedef const_iterator1 iterator1;
0406     class const_iterator2;
0407     typedef const_iterator2 iterator2;
0408 #endif
0409     typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
0410     typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
0411
0412     // Element lookup
0413     BOOST_UBLAS_INLINE
0414     const_iterator1 find1 (int rank, size_type i, size_type j) const {
0415         const_subiterator1_type it1 (e1_.find (i));
0416         const_subiterator1_type it1_end (e1_.find (size1 ()));
0417         const_subiterator2_type it2 (e2_.find (j));
0418         const_subiterator2_type it2_end (e2_.find (size2 ()));
0419         if (it2 == it2_end || (rank == 1 && (it2.index () != j || *it2 == typename E2::value_type/*zero*/()))) {
0420             it1 = it1_end;
0421             it2 = it2_end;
0422         }
0423 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
0424         return const_iterator1 (*this, it1.index (), it2.index ());
0425 #else
0426 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
0427         return const_iterator1 (*this, it1, it2, it2 != it2_end ? *it2 : typename E2::value_type/*zero*/());
0428 #else
0429         return const_iterator1 (*this, it1, it2);
0430 #endif
0431 #endif
0432     }
0433     BOOST_UBLAS_INLINE
0434     const_iterator2 find2 (int rank, size_type i, size_type j) const {
0435         const_subiterator2_type it2 (e2_.find (j));
0436         const_subiterator2_type it2_end (e2_.find (size2 ()));
0437         const_subiterator1_type it1 (e1_.find (i));
0438         const_subiterator1_type it1_end (e1_.find (size1 ()));
0439         if (it1 == it1_end || (rank == 1 && (it1.index () != i || *it1 == typename E1::value_type/*zero*/()))) {
0440             it2 = it2_end;
0441             it1 = it1_end;
0442         }
0443 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
0444         return const_iterator2 (*this, it1.index (), it2.index ());
0445 #else
0446 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
0447         return const_iterator2 (*this, it1, it2, it1 != it1_end ? *it1 : typename E1::value_type/*zero*/());
0448 #else
0449         return const_iterator2 (*this, it1, it2);
0450 #endif
0451 #endif
0452     }
0453
0454     // Iterators enhance the iterators of the referenced expressions
0455     // with the binary functor.
0456
0457 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
0458     class const_iterator1:
0459         public container_const_reference<vector_matrix_binary>,
0460         public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
0461         typename E2::const_iterator::iterator_category>::iterator_category>::template
0462         iterator_base<const_iterator1, value_type>::type {
0463     public:
0464         typedef typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
0465         typename E2::const_iterator::iterator_category>::iterator_category iterator_category;
0466         typedef typename vector_matrix_binary::difference_type difference_type;
0467         typedef typename vector_matrix_binary::value_type value_type;
0468         typedef typename vector_matrix_binary::const_reference reference;
0469         typedef typename vector_matrix_binary::const_pointer pointer;
0470
0471         typedef const_iterator2 dual_iterator_type;
0472         typedef const_reverse_iterator2 dual_reverse_iterator_type;
0473
0474         // Construction and destruction
0475 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
0476         BOOST_UBLAS_INLINE
0477         const_iterator1 ():
0478           container_const_reference<self_type> (), it1_ (), it2_ (), t2_ () {}
0479         BOOST_UBLAS_INLINE
0480         const_iterator1 (const self_type &vmb, const const_subiterator1_type &it1, const const_subiterator2_type &it2, value_type t2):
0481           container_const_reference<self_type> (vmb), it1_ (it1), it2_ (it2), t2_ (t2) {}
0482 #else
0483         BOOST_UBLAS_INLINE
0484         const_iterator1 ():
0485           container_const_reference<self_type> (), it1_ (), it2_ () {}
0486         BOOST_UBLAS_INLINE
0487         const_iterator1 (const self_type &vmb, const const_subiterator1_type &it1, const const_subiterator2_type &it2):
0488           container_const_reference<self_type> (vmb), it1_ (it1), it2_ (it2) {}
0489 #endif
0490
0491         // Arithmetic
0492         BOOST_UBLAS_INLINE
0493         const_iterator1 &operator ++ () {
0494             ++ it1_;
0495             return *this;
0496         }
0497         BOOST_UBLAS_INLINE
0498         const_iterator1 &operator -- () {
0499             -- it1_;
0500             return *this;
0501         }
0502         BOOST_UBLAS_INLINE
0503         const_iterator1 &operator += (difference_type n) {
0504             it1_ += n;
0505             return *this;
0506         }
0507         BOOST_UBLAS_INLINE
0508         const_iterator1 &operator -= (difference_type n) {
0509             it1_ -= n;
0510             return *this;
0511         }
0512         BOOST_UBLAS_INLINE
0513         difference_type operator - (const const_iterator1 &it) const {
0514             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
0515             BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
0516             return it1_ - it.it1_;
0517         }
0518
0519         // Dereference
0520         BOOST_UBLAS_INLINE
0521         const_reference operator * () const {
0522 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
0523             return functor_type::apply (*it1_, t2_);
0524 #else
0525             return functor_type::apply (*it1_, *it2_);
0526 #endif
0527         }
0528         BOOST_UBLAS_INLINE
0529         const_reference operator [] (difference_type n) const {
0530             return *(*this + n);
0531         }
0532
0533 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
0534         BOOST_UBLAS_INLINE
0535 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
0536         typename self_type::
0537 #endif
0538         const_iterator2 begin () const {
0539             return (*this) ().find2 (1, index1 (), 0);
0540         }
0541         BOOST_UBLAS_INLINE
0542 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
0543         typename self_type::
0544 #endif
0545         const_iterator2 cbegin () const {
0546             return begin ();
0547         }
0548         BOOST_UBLAS_INLINE
0549 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
0550         typename self_type::
0551 #endif
0552         const_iterator2 end () const {
0553             return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
0554         }
0555         BOOST_UBLAS_INLINE
0556 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
0557         typename self_type::
0558 #endif
0559         const_iterator2 cend () const {
0560             return end ();
0561         }
0562         BOOST_UBLAS_INLINE
0563 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
0564         typename self_type::
0565 #endif
0566         const_reverse_iterator2 rbegin () const {
0567             return const_reverse_iterator2 (end ());
0568         }
0569         BOOST_UBLAS_INLINE
0570 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
0571         typename self_type::
0572 #endif
0573         const_reverse_iterator2 crbegin () const {
0574             return rbegin ();
0575         }
0576         BOOST_UBLAS_INLINE
0577 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
0578         typename self_type::
0579 #endif
0580         const_reverse_iterator2 rend () const {
0581             return const_reverse_iterator2 (begin ());
0582         }
0583         BOOST_UBLAS_INLINE
0584 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
0585         typename self_type::
0586 #endif
0587         const_reverse_iterator2 crend () const {
0588             return rend ();
0589         }
0590 #endif
0591
0592         // Indices
0593         BOOST_UBLAS_INLINE
0594         size_type index1 () const {
0595             return it1_.index ();
0596         }
0597         BOOST_UBLAS_INLINE
0598         size_type  index2 () const {
0599             return it2_.index ();
0600         }
0601
0602         // Assignment
0603         BOOST_UBLAS_INLINE
0604         const_iterator1 &operator = (const const_iterator1 &it) {
0605             container_const_reference<self_type>::assign (&it ());
0606             it1_ = it.it1_;
0607             it2_ = it.it2_;
0608 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
0609             t2_ = it.t2_;
0610 #endif
0611             return *this;
0612         }
0613
0614         // Comparison
0615         BOOST_UBLAS_INLINE
0616         bool operator == (const const_iterator1 &it) const {
0617             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
0618             BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
0619             return it1_ == it.it1_;
0620         }
0621         BOOST_UBLAS_INLINE
0622         bool operator < (const const_iterator1 &it) const {
0623             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
0624             BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
0625             return it1_ < it.it1_;
0626         }
0627
0628     private:
0629 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
0630         const_subiterator1_type it1_;
0631         // Mutable due to assignment
0632         /* const */ const_subiterator2_type it2_;
0633         value_type t2_;
0634 #else
0635         const_subiterator1_type it1_;
0636         const_subiterator2_type it2_;
0637 #endif
0638     };
0639 #endif
0640
0641     BOOST_UBLAS_INLINE
0642     const_iterator1 begin1 () const {
0643         return find1 (0, 0, 0);
0644     }
0645     BOOST_UBLAS_INLINE
0646     const_iterator1 cbegin1 () const {
0647         return begin1 ();
0648     }
0649     BOOST_UBLAS_INLINE
0650     const_iterator1 end1 () const {
0651         return find1 (0, size1 (), 0);
0652     }
0653     BOOST_UBLAS_INLINE
0654     const_iterator1 cend1 () const {
0655         return end1 ();
0656     }
0657
0658 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
0659     class const_iterator2:
0660         public container_const_reference<vector_matrix_binary>,
0661         public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
0662         typename E2::const_iterator::iterator_category>::iterator_category>::template
0663         iterator_base<const_iterator2, value_type>::type {
0664     public:
0665         typedef typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
0666         typename E2::const_iterator::iterator_category>::iterator_category iterator_category;
0667         typedef typename vector_matrix_binary::difference_type difference_type;
0668         typedef typename vector_matrix_binary::value_type value_type;
0669         typedef typename vector_matrix_binary::const_reference reference;
0670         typedef typename vector_matrix_binary::const_pointer pointer;
0671
0672         typedef const_iterator1 dual_iterator_type;
0673         typedef const_reverse_iterator1 dual_reverse_iterator_type;
0674
0675         // Construction and destruction
0676 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
0677         BOOST_UBLAS_INLINE
0678         const_iterator2 ():
0679           container_const_reference<self_type> (), it1_ (), it2_ (), t1_ () {}
0680         BOOST_UBLAS_INLINE
0681         const_iterator2 (const self_type &vmb, const const_subiterator1_type &it1, const const_subiterator2_type &it2, value_type t1):
0682           container_const_reference<self_type> (vmb), it1_ (it1), it2_ (it2), t1_ (t1) {}
0683 #else
0684         BOOST_UBLAS_INLINE
0685         const_iterator2 ():
0686           container_const_reference<self_type> (), it1_ (), it2_ () {}
0687         BOOST_UBLAS_INLINE
0688         const_iterator2 (const self_type &vmb, const const_subiterator1_type &it1, const const_subiterator2_type &it2):
0689           container_const_reference<self_type> (vmb), it1_ (it1), it2_ (it2) {}
0690 #endif
0691
0692         // Arithmetic
0693         BOOST_UBLAS_INLINE
0694         const_iterator2 &operator ++ () {
0695             ++ it2_;
0696             return *this;
0697         }
0698         BOOST_UBLAS_INLINE
0699         const_iterator2 &operator -- () {
0700             -- it2_;
0701             return *this;
0702         }
0703         BOOST_UBLAS_INLINE
0704         const_iterator2 &operator += (difference_type n) {
0705             it2_ += n;
0706             return *this;
0707         }
0708         BOOST_UBLAS_INLINE
0709         const_iterator2 &operator -= (difference_type n) {
0710             it2_ -= n;
0711             return *this;
0712         }
0713         BOOST_UBLAS_INLINE
0714         difference_type operator - (const const_iterator2 &it) const {
0715             BOOST_UBLAS_CHECK ((*this) ().same_closure(it ()), external_logic ());
0716             BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
0717             return it2_ - it.it2_;
0718         }
0719
0720         // Dereference
0721         BOOST_UBLAS_INLINE
0722         const_reference operator * () const {
0723 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
0724             return functor_type::apply (t1_, *it2_);
0725 #else
0726             return functor_type::apply (*it1_, *it2_);
0727 #endif
0728         }
0729         BOOST_UBLAS_INLINE
0730         const_reference operator [] (difference_type n) const {
0731             return *(*this + n);
0732         }
0733
0734 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
0735         BOOST_UBLAS_INLINE
0736 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
0737         typename self_type::
0738 #endif
0739         const_iterator1 begin () const {
0740             return (*this) ().find1 (1, 0, index2 ());
0741         }
0742         BOOST_UBLAS_INLINE
0743 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
0744         typename self_type::
0745 #endif
0746         const_iterator1 cbegin () const {
0747             return begin ();
0748         }
0749         BOOST_UBLAS_INLINE
0750 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
0751         typename self_type::
0752 #endif
0753         const_iterator1 end () const {
0754             return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
0755         }
0756         BOOST_UBLAS_INLINE
0757 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
0758         typename self_type::
0759 #endif
0760         const_iterator1 cend () const {
0761             return end ();
0762         }
0763         BOOST_UBLAS_INLINE
0764 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
0765         typename self_type::
0766 #endif
0767         const_reverse_iterator1 rbegin () const {
0768             return const_reverse_iterator1 (end ());
0769         }
0770         BOOST_UBLAS_INLINE
0771 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
0772         typename self_type::
0773 #endif
0774         const_reverse_iterator1 crbegin () const {
0775             return rbegin ();
0776         }
0777         BOOST_UBLAS_INLINE
0778 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
0779         typename self_type::
0780 #endif
0781         const_reverse_iterator1 rend () const {
0782             return const_reverse_iterator1 (begin ());
0783         }
0784         BOOST_UBLAS_INLINE
0785 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
0786         typename self_type::
0787 #endif
0788         const_reverse_iterator1 crend () const {
0789             return rend ();
0790         }
0791 #endif
0792
0793         // Indices
0794         BOOST_UBLAS_INLINE
0795         size_type index1 () const {
0796             return it1_.index ();
0797         }
0798         BOOST_UBLAS_INLINE
0799         size_type  index2 () const {
0800             return it2_.index ();
0801         }
0802
0803         // Assignment
0804         BOOST_UBLAS_INLINE
0805         const_iterator2 &operator = (const const_iterator2 &it) {
0806             container_const_reference<self_type>::assign (&it ());
0807             it1_ = it.it1_;
0808             it2_ = it.it2_;
0809 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
0810             t1_ = it.t1_;
0811 #endif
0812             return *this;
0813         }
0814
0815         // Comparison
0816         BOOST_UBLAS_INLINE
0817         bool operator == (const const_iterator2 &it) const {
0818             BOOST_UBLAS_CHECK ((*this) ().same_closure( it ()), external_logic ());
0819             BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
0820             return it2_ == it.it2_;
0821         }
0822         BOOST_UBLAS_INLINE
0823         bool operator < (const const_iterator2 &it) const {
0824             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
0825             BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
0826             return it2_ < it.it2_;
0827         }
0828
0829     private:
0830 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
0831         // Mutable due to assignment
0832         /* const */ const_subiterator1_type it1_;
0833         const_subiterator2_type it2_;
0834         value_type t1_;
0835 #else
0836         const_subiterator1_type it1_;
0837         const_subiterator2_type it2_;
0838 #endif
0839     };
0840 #endif
0841
0842     BOOST_UBLAS_INLINE
0843     const_iterator2 begin2 () const {
0844         return find2 (0, 0, 0);
0845     }
0846     BOOST_UBLAS_INLINE
0847     const_iterator2 cbegin2 () const {
0848         return begin2 ();
0849     }
0850     BOOST_UBLAS_INLINE
0851     const_iterator2 end2 () const {
0852         return find2 (0, 0, size2 ());
0853     }
0854     BOOST_UBLAS_INLINE
0855     const_iterator2 cend2 () const {
0856         return end2 ();
0857     }
0858
0859     // Reverse iterators
0860
0861     BOOST_UBLAS_INLINE
0862     const_reverse_iterator1 rbegin1 () const {
0863         return const_reverse_iterator1 (end1 ());
0864     }
0865     BOOST_UBLAS_INLINE
0866     const_reverse_iterator1 crbegin1 () const {
0867         return rbegin1 ();
0868     }
0869     BOOST_UBLAS_INLINE
0870     const_reverse_iterator1 rend1 () const {
0871         return const_reverse_iterator1 (begin1 ());
0872     }
0873     BOOST_UBLAS_INLINE
0874     const_reverse_iterator1 crend1 () const {
0875         return rend1 ();
0876     }
0877     BOOST_UBLAS_INLINE
0878     const_reverse_iterator2 rbegin2 () const {
0879         return const_reverse_iterator2 (end2 ());
0880     }
0881     BOOST_UBLAS_INLINE
0882     const_reverse_iterator2 crbegin2 () const {
0883         return rbegin2 ();
0884     }
0885     BOOST_UBLAS_INLINE
0886     const_reverse_iterator2 rend2 () const {
0887         return const_reverse_iterator2 (begin2 ());
0888     }
0889     BOOST_UBLAS_INLINE
0890     const_reverse_iterator2 crend2 () const {
0891         return rend2 ();
0892     }
0893
0894 private:
0895     expression1_closure_type e1_;
0896     expression2_closure_type e2_;
0897 };
0898
0899 template<class E1, class E2, class F>
0900 struct vector_matrix_binary_traits {
0901     typedef vector_matrix_binary<E1, E2, F> expression_type;
0902 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
0903     typedef expression_type result_type;
0904 #else
0905     // ISSUE matrix is arbitary temporary type
0906     typedef matrix<typename F::value_type> result_type;
0907 #endif
0908 };
0909
0910 // (outer_prod (v1, v2)) [i] [j] = v1 [i] * v2 [j]
0911 template<class E1, class E2>
0912 BOOST_UBLAS_INLINE
0913 typename vector_matrix_binary_traits<E1, E2, scalar_multiplies<typename E1::value_type, typename E2::value_type> >::result_type
0914 outer_prod (const vector_expression<E1> &e1,
0915             const vector_expression<E2> &e2) {
0916     BOOST_STATIC_ASSERT (E1::complexity == 0 && E2::complexity == 0);
0917     typedef typename vector_matrix_binary_traits<E1, E2, scalar_multiplies<typename E1::value_type, typename E2::value_type> >::expression_type expression_type;
0918     return expression_type (e1 (), e2 ());
0919 }
0920
0921 template<class E, class F>
0922 class matrix_unary1:
0923     public matrix_expression<matrix_unary1<E, F> > {
0924
0925     typedef E expression_type;
0926     typedef F functor_type;
0927 public:
0928     typedef typename E::const_closure_type expression_closure_type;
0929 private:
0930     typedef matrix_unary1<E, F> self_type;
0931 public:
0932 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
0933     using matrix_expression<self_type>::operator ();
0934 #endif
0935     typedef typename E::size_type size_type;
0936     typedef typename E::difference_type difference_type;
0937     typedef typename F::result_type value_type;
0938     typedef value_type const_reference;
0939     typedef const_reference reference;
0940     typedef const self_type const_closure_type;
0941     typedef const_closure_type closure_type;
0942     typedef typename E::orientation_category orientation_category;
0943     typedef unknown_storage_tag storage_category;
0944
0945     // Construction and destruction
0946     BOOST_UBLAS_INLINE
0947     explicit matrix_unary1 (const expression_type &e):
0948       e_ (e) {}
0949
0950     // Accessors
0951     BOOST_UBLAS_INLINE
0952     size_type size1 () const {
0953         return e_.size1 ();
0954     }
0955     BOOST_UBLAS_INLINE
0956     size_type size2 () const {
0957         return e_.size2 ();
0958     }
0959
0960 public:
0961     // Expression accessors
0962     BOOST_UBLAS_INLINE
0963     const expression_closure_type &expression () const {
0964         return e_;
0965     }
0966
0967 public:
0968     // Element access
0969     BOOST_UBLAS_INLINE
0970     const_reference operator () (size_type i, size_type j) const {
0971         return functor_type::apply (e_ (i, j));
0972     }
0973
0974     // Element access
0975     BOOST_UBLAS_INLINE
0976     const_reference operator () (size_type i) const {
0977         return functor_type::apply (e_ (i));
0978     }
0979
0980
0981     // Closure comparison
0982     BOOST_UBLAS_INLINE
0983     bool same_closure (const matrix_unary1 &mu1) const {
0984         return (*this).expression ().same_closure (mu1.expression ());
0985     }
0986
0987     // Iterator types
0988 private:
0989     typedef typename E::const_iterator1 const_subiterator1_type;
0990     typedef typename E::const_iterator2 const_subiterator2_type;
0991     typedef const value_type *const_pointer;
0992
0993 public:
0994 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
0995     typedef indexed_const_iterator1<const_closure_type, typename const_subiterator1_type::iterator_category> const_iterator1;
0996     typedef const_iterator1 iterator1;
0997     typedef indexed_const_iterator2<const_closure_type, typename const_subiterator2_type::iterator_category> const_iterator2;
0998     typedef const_iterator2 iterator2;
0999 #else
1000     class const_iterator1;
1001     typedef const_iterator1 iterator1;
1002     class const_iterator2;
1003     typedef const_iterator2 iterator2;
1004 #endif
1005     typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
1006     typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
1007
1008     // Element lookup
1009     BOOST_UBLAS_INLINE
1010     const_iterator1 find1 (int rank, size_type i, size_type j) const {
1011         const_subiterator1_type it1 (e_.find1 (rank, i, j));
1012 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1013         return const_iterator1 (*this, it1.index1 (), it1.index2 ());
1014 #else
1015         return const_iterator1 (*this, it1);
1016 #endif
1017     }
1018     BOOST_UBLAS_INLINE
1019     const_iterator2 find2 (int rank, size_type i, size_type j) const {
1020         const_subiterator2_type it2 (e_.find2 (rank, i, j));
1021 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1022         return const_iterator2 (*this, it2.index1 (), it2.index2 ());
1023 #else
1024         return const_iterator2 (*this, it2);
1025 #endif
1026     }
1027
1028     // Iterators enhance the iterators of the referenced expression
1029     // with the unary functor.
1030
1031 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1032     class const_iterator1:
1033         public container_const_reference<matrix_unary1>,
1034         public iterator_base_traits<typename E::const_iterator1::iterator_category>::template
1035         iterator_base<const_iterator1, value_type>::type {
1036     public:
1037         typedef typename E::const_iterator1::iterator_category iterator_category;
1038         typedef typename matrix_unary1::difference_type difference_type;
1039         typedef typename matrix_unary1::value_type value_type;
1040         typedef typename matrix_unary1::const_reference reference;
1041         typedef typename matrix_unary1::const_pointer pointer;
1042
1043         typedef const_iterator2 dual_iterator_type;
1044         typedef const_reverse_iterator2 dual_reverse_iterator_type;
1045
1046         // Construction and destruction
1047         BOOST_UBLAS_INLINE
1048         const_iterator1 ():
1049           container_const_reference<self_type> (), it_ () {}
1050         BOOST_UBLAS_INLINE
1051         const_iterator1 (const self_type &mu, const const_subiterator1_type &it):
1052           container_const_reference<self_type> (mu), it_ (it) {}
1053
1054         // Arithmetic
1055         BOOST_UBLAS_INLINE
1056         const_iterator1 &operator ++ () {
1057             ++ it_;
1058             return *this;
1059         }
1060         BOOST_UBLAS_INLINE
1061         const_iterator1 &operator -- () {
1062             -- it_;
1063             return *this;
1064         }
1065         BOOST_UBLAS_INLINE
1066         const_iterator1 &operator += (difference_type n) {
1067             it_ += n;
1068             return *this;
1069         }
1070         BOOST_UBLAS_INLINE
1071         const_iterator1 &operator -= (difference_type n) {
1072             it_ -= n;
1073             return *this;
1074         }
1075         BOOST_UBLAS_INLINE
1076         difference_type operator - (const const_iterator1 &it) const {
1077             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1078             return it_ - it.it_;
1079         }
1080
1081         // Dereference
1082         BOOST_UBLAS_INLINE
1083         const_reference operator * () const {
1084             return functor_type::apply (*it_);
1085         }
1086         BOOST_UBLAS_INLINE
1087         const_reference operator [] (difference_type n) const {
1088             return *(*this + n);
1089         }
1090
1091 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1092         BOOST_UBLAS_INLINE
1093 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1094         typename self_type::
1095 #endif
1096         const_iterator2 begin () const {
1097             return (*this) ().find2 (1, index1 (), 0);
1098         }
1099         BOOST_UBLAS_INLINE
1100 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1101         typename self_type::
1102 #endif
1103         const_iterator2 cbegin () const {
1104             return begin ();
1105         }
1106         BOOST_UBLAS_INLINE
1107 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1108         typename self_type::
1109 #endif
1110         const_iterator2 end () const {
1111             return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
1112         }
1113         BOOST_UBLAS_INLINE
1114 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1115         typename self_type::
1116 #endif
1117         const_iterator2 cend () const {
1118             return end ();
1119         }
1120         BOOST_UBLAS_INLINE
1121 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1122         typename self_type::
1123 #endif
1124         const_reverse_iterator2 rbegin () const {
1125             return const_reverse_iterator2 (end ());
1126         }
1127         BOOST_UBLAS_INLINE
1128 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1129         typename self_type::
1130 #endif
1131         const_reverse_iterator2 crbegin () const {
1132             return rbegin ();
1133         }
1134         BOOST_UBLAS_INLINE
1135 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1136         typename self_type::
1137 #endif
1138         const_reverse_iterator2 rend () const {
1139             return const_reverse_iterator2 (begin ());
1140         }
1141         BOOST_UBLAS_INLINE
1142 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1143         typename self_type::
1144 #endif
1145         const_reverse_iterator2 crend () const {
1146             return rend ();
1147         }
1148 #endif
1149
1150         // Indices
1151         BOOST_UBLAS_INLINE
1152         size_type index1 () const {
1153             return it_.index1 ();
1154         }
1155         BOOST_UBLAS_INLINE
1156         size_type index2 () const {
1157             return it_.index2 ();
1158         }
1159
1160         // Assignment
1161         BOOST_UBLAS_INLINE
1162         const_iterator1 &operator = (const const_iterator1 &it) {
1163             container_const_reference<self_type>::assign (&it ());
1164             it_ = it.it_;
1165             return *this;
1166         }
1167
1168         // Comparison
1169         BOOST_UBLAS_INLINE
1170         bool operator == (const const_iterator1 &it) const {
1171             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1172             return it_ == it.it_;
1173         }
1174         BOOST_UBLAS_INLINE
1175         bool operator < (const const_iterator1 &it) const {
1176             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1177             return it_ < it.it_;
1178         }
1179
1180     private:
1181         const_subiterator1_type it_;
1182     };
1183 #endif
1184
1185     BOOST_UBLAS_INLINE
1186     const_iterator1 begin1 () const {
1187         return find1 (0, 0, 0);
1188     }
1189     BOOST_UBLAS_INLINE
1190     const_iterator1 cbegin1 () const {
1191         return begin1 ();
1192     }
1193     BOOST_UBLAS_INLINE
1194     const_iterator1 end1 () const {
1195         return find1 (0, size1 (), 0);
1196     }
1197     BOOST_UBLAS_INLINE
1198     const_iterator1 cend1 () const {
1199         return end1 ();
1200     }
1201
1202 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1203     class const_iterator2:
1204         public container_const_reference<matrix_unary1>,
1205         public iterator_base_traits<typename E::const_iterator2::iterator_category>::template
1206         iterator_base<const_iterator2, value_type>::type {
1207     public:
1208         typedef typename E::const_iterator2::iterator_category iterator_category;
1209         typedef typename matrix_unary1::difference_type difference_type;
1210         typedef typename matrix_unary1::value_type value_type;
1211         typedef typename matrix_unary1::const_reference reference;
1212         typedef typename matrix_unary1::const_pointer pointer;
1213
1214         typedef const_iterator1 dual_iterator_type;
1215         typedef const_reverse_iterator1 dual_reverse_iterator_type;
1216
1217         // Construction and destruction
1218         BOOST_UBLAS_INLINE
1219         const_iterator2 ():
1220           container_const_reference<self_type> (), it_ () {}
1221         BOOST_UBLAS_INLINE
1222         const_iterator2 (const self_type &mu, const const_subiterator2_type &it):
1223           container_const_reference<self_type> (mu), it_ (it) {}
1224
1225         // Arithmetic
1226         BOOST_UBLAS_INLINE
1227         const_iterator2 &operator ++ () {
1228             ++ it_;
1229             return *this;
1230         }
1231         BOOST_UBLAS_INLINE
1232         const_iterator2 &operator -- () {
1233             -- it_;
1234             return *this;
1235         }
1236         BOOST_UBLAS_INLINE
1237         const_iterator2 &operator += (difference_type n) {
1238             it_ += n;
1239             return *this;
1240         }
1241         BOOST_UBLAS_INLINE
1242         const_iterator2 &operator -= (difference_type n) {
1243             it_ -= n;
1244             return *this;
1245         }
1246         BOOST_UBLAS_INLINE
1247         difference_type operator - (const const_iterator2 &it) const {
1248             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1249             return it_ - it.it_;
1250         }
1251
1252         // Dereference
1253         BOOST_UBLAS_INLINE
1254         const_reference operator * () const {
1255             return functor_type::apply (*it_);
1256         }
1257         BOOST_UBLAS_INLINE
1258         const_reference operator [] (difference_type n) const {
1259             return *(*this + n);
1260         }
1261
1262 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1263         BOOST_UBLAS_INLINE
1264 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1265         typename self_type::
1266 #endif
1267         const_iterator1 begin () const {
1268             return (*this) ().find1 (1, 0, index2 ());
1269         }
1270         BOOST_UBLAS_INLINE
1271 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1272         typename self_type::
1273 #endif
1274         const_iterator1 cbegin () const {
1275             return begin ();
1276         }
1277         BOOST_UBLAS_INLINE
1278 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1279         typename self_type::
1280 #endif
1281         const_iterator1 end () const {
1282             return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
1283         }
1284         BOOST_UBLAS_INLINE
1285 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1286         typename self_type::
1287 #endif
1288         const_iterator1 cend () const {
1289             return end ();
1290         }
1291         BOOST_UBLAS_INLINE
1292 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1293         typename self_type::
1294 #endif
1295         const_reverse_iterator1 rbegin () const {
1296             return const_reverse_iterator1 (end ());
1297         }
1298         BOOST_UBLAS_INLINE
1299 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1300         typename self_type::
1301 #endif
1302         const_reverse_iterator1 crbegin () const {
1303             return rbegin ();
1304         }
1305         BOOST_UBLAS_INLINE
1306 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1307         typename self_type::
1308 #endif
1309         const_reverse_iterator1 rend () const {
1310             return const_reverse_iterator1 (begin ());
1311         }
1312         BOOST_UBLAS_INLINE
1313 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1314         typename self_type::
1315 #endif
1316         const_reverse_iterator1 crend () const {
1317             return rend ();
1318         }
1319 #endif
1320
1321         // Indices
1322         BOOST_UBLAS_INLINE
1323         size_type index1 () const {
1324             return it_.index1 ();
1325         }
1326         BOOST_UBLAS_INLINE
1327         size_type index2 () const {
1328             return it_.index2 ();
1329         }
1330
1331         // Assignment
1332         BOOST_UBLAS_INLINE
1333         const_iterator2 &operator = (const const_iterator2 &it) {
1334             container_const_reference<self_type>::assign (&it ());
1335             it_ = it.it_;
1336             return *this;
1337         }
1338
1339         // Comparison
1340         BOOST_UBLAS_INLINE
1341         bool operator == (const const_iterator2 &it) const {
1342             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1343             return it_ == it.it_;
1344         }
1345         BOOST_UBLAS_INLINE
1346         bool operator < (const const_iterator2 &it) const {
1347             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1348             return it_ < it.it_;
1349         }
1350
1351     private:
1352         const_subiterator2_type it_;
1353     };
1354 #endif
1355
1356     BOOST_UBLAS_INLINE
1357     const_iterator2 begin2 () const {
1358         return find2 (0, 0, 0);
1359     }
1360     BOOST_UBLAS_INLINE
1361     const_iterator2 cbegin2 () const {
1362         return begin2 ();
1363     }
1364     BOOST_UBLAS_INLINE
1365     const_iterator2 end2 () const {
1366         return find2 (0, 0, size2 ());
1367     }
1368     BOOST_UBLAS_INLINE
1369     const_iterator2 cend2 () const {
1370         return end2 ();
1371     }
1372
1373     // Reverse iterators
1374
1375     BOOST_UBLAS_INLINE
1376     const_reverse_iterator1 rbegin1 () const {
1377         return const_reverse_iterator1 (end1 ());
1378     }
1379     BOOST_UBLAS_INLINE
1380     const_reverse_iterator1 crbegin1 () const {
1381         return rbegin1 ();
1382     }
1383     BOOST_UBLAS_INLINE
1384     const_reverse_iterator1 rend1 () const {
1385         return const_reverse_iterator1 (begin1 ());
1386     }
1387     BOOST_UBLAS_INLINE
1388     const_reverse_iterator1 crend1 () const {
1389         return rend1 ();
1390     }
1391
1392     BOOST_UBLAS_INLINE
1393     const_reverse_iterator2 rbegin2 () const {
1394         return const_reverse_iterator2 (end2 ());
1395     }
1396     BOOST_UBLAS_INLINE
1397     const_reverse_iterator2 crbegin2 () const {
1398         return rbegin2 ();
1399     }
1400     BOOST_UBLAS_INLINE
1401     const_reverse_iterator2 rend2 () const {
1402         return const_reverse_iterator2 (begin2 ());
1403     }
1404     BOOST_UBLAS_INLINE
1405     const_reverse_iterator2 crend2 () const {
1406         return rend2 ();
1407     }
1408
1409 private:
1410     expression_closure_type e_;
1411 };
1412
1413 template<class E, class F>
1414 struct matrix_unary1_traits {
1415     typedef matrix_unary1<E, F> expression_type;
1416 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
1417     typedef expression_type result_type;
1418 #else
1419     typedef typename E::matrix_temporary_type result_type;
1420 #endif
1421 };
1422
1423 // (- m) [i] [j] = - m [i] [j]
1424 template<class E>
1425 BOOST_UBLAS_INLINE
1426 typename matrix_unary1_traits<E, scalar_negate<typename E::value_type> >::result_type
1427 operator - (const matrix_expression<E> &e) {
1428     typedef typename matrix_unary1_traits<E, scalar_negate<typename E::value_type> >::expression_type expression_type;
1429     return expression_type (e ());
1430 }
1431
1432 // (conj m) [i] [j] = conj (m [i] [j])
1433 template<class E>
1434 BOOST_UBLAS_INLINE
1435 typename matrix_unary1_traits<E, scalar_conj<typename E::value_type> >::result_type
1436 conj (const matrix_expression<E> &e) {
1437     typedef typename matrix_unary1_traits<E, scalar_conj<typename E::value_type> >::expression_type expression_type;
1438     return expression_type (e ());
1439 }
1440
1441 // (real m) [i] [j] = real (m [i] [j])
1442 template<class E>
1443 BOOST_UBLAS_INLINE
1444 typename matrix_unary1_traits<E, scalar_real<typename E::value_type> >::result_type
1445 real (const matrix_expression<E> &e) {
1446     typedef typename matrix_unary1_traits<E, scalar_real<typename E::value_type> >::expression_type expression_type;
1447     return expression_type (e ());
1448 }
1449
1450 // (imag m) [i] [j] = imag (m [i] [j])
1451 template<class E>
1452 BOOST_UBLAS_INLINE
1453 typename matrix_unary1_traits<E, scalar_imag<typename E::value_type> >::result_type
1454 imag (const matrix_expression<E> &e) {
1455     typedef typename matrix_unary1_traits<E, scalar_imag<typename E::value_type> >::expression_type expression_type;
1456     return expression_type (e ());
1457 }
1458
1459 template<class E, class F>
1460 class matrix_unary2:
1461     public matrix_expression<matrix_unary2<E, F> > {
1462
1463     typedef typename boost::mpl::if_<boost::is_same<F, scalar_identity<typename E::value_type> >,
1464     E,
1465     const E>::type expression_type;
1466     typedef F functor_type;
1467 public:
1468     typedef typename boost::mpl::if_<boost::is_const<expression_type>,
1469     typename E::const_closure_type,
1470     typename E::closure_type>::type expression_closure_type;
1471 private:
1472     typedef matrix_unary2<E, F> self_type;
1473 public:
1474 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
1475     using matrix_expression<self_type>::operator ();
1476 #endif
1477     typedef typename E::size_type size_type;
1478     typedef typename E::difference_type difference_type;
1479     typedef typename F::result_type value_type;
1480     typedef value_type const_reference;
1481     typedef typename boost::mpl::if_<boost::is_same<F, scalar_identity<value_type> >,
1482     typename E::reference,
1483     value_type>::type reference;
1484
1485     typedef const self_type const_closure_type;
1486     typedef self_type closure_type;
1487     typedef typename boost::mpl::if_<boost::is_same<typename E::orientation_category,
1488     row_major_tag>,
1489     column_major_tag,
1490     typename boost::mpl::if_<boost::is_same<typename E::orientation_category,
1491     column_major_tag>,
1492     row_major_tag,
1493     typename E::orientation_category>::type>::type orientation_category;
1494     typedef typename E::storage_category storage_category;
1495
1496     // Construction and destruction
1497     BOOST_UBLAS_INLINE
1498     // matrix_unary2 may be used as mutable expression -
1499     // this is the only non const expression constructor
1500     explicit matrix_unary2 (expression_type &e):
1501       e_ (e) {}
1502
1503     // Accessors
1504     BOOST_UBLAS_INLINE
1505     size_type size1 () const {
1506         return e_.size2 ();
1507     }
1508     BOOST_UBLAS_INLINE
1509     size_type size2 () const {
1510         return e_.size1 ();
1511     }
1512
1513 public:
1514     // Expression accessors
1515     BOOST_UBLAS_INLINE
1516     const expression_closure_type &expression () const {
1517         return e_;
1518     }
1519
1520 public:
1521     // Element access
1522     BOOST_UBLAS_INLINE
1523     const_reference operator () (size_type i, size_type j) const {
1524         return functor_type::apply (e_ (j, i));
1525     }
1526     BOOST_UBLAS_INLINE
1527     reference operator () (size_type i, size_type j) {
1528         BOOST_STATIC_ASSERT ((boost::is_same<functor_type, scalar_identity<value_type > >::value));
1529         return e_ (j, i);
1530     }
1531
1532     // Element access
1533     BOOST_UBLAS_INLINE
1534     const_reference operator () (size_type i) const {
1535         return functor_type::apply (e_ (i));
1536     }
1537     BOOST_UBLAS_INLINE
1538     reference operator () (size_type i) {
1539         BOOST_STATIC_ASSERT ((boost::is_same<functor_type, scalar_identity<value_type > >::value));
1540         return e_ (i);
1541     }
1542
1543     // Closure comparison
1544     BOOST_UBLAS_INLINE
1545     bool same_closure (const matrix_unary2 &mu2) const {
1546         return (*this).expression ().same_closure (mu2.expression ());
1547     }
1548
1549     // Iterator types
1550 private:
1551     typedef typename E::const_iterator1 const_subiterator2_type;
1552     typedef typename E::const_iterator2 const_subiterator1_type;
1553     typedef const value_type *const_pointer;
1554
1555 public:
1556 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1557     typedef indexed_const_iterator1<const_closure_type, typename const_subiterator1_type::iterator_category> const_iterator1;
1558     typedef const_iterator1 iterator1;
1559     typedef indexed_const_iterator2<const_closure_type, typename const_subiterator2_type::iterator_category> const_iterator2;
1560     typedef const_iterator2 iterator2;
1561 #else
1562     class const_iterator1;
1563     typedef const_iterator1 iterator1;
1564     class const_iterator2;
1565     typedef const_iterator2 iterator2;
1566 #endif
1567     typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
1568     typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
1569
1570     // Element lookup
1571     BOOST_UBLAS_INLINE
1572     const_iterator1 find1 (int rank, size_type i, size_type j) const {
1573         const_subiterator1_type it1 (e_.find2 (rank, j, i));
1574 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1575         return const_iterator1 (*this, it1.index2 (), it1.index1 ());
1576 #else
1577         return const_iterator1 (*this, it1);
1578 #endif
1579     }
1580     BOOST_UBLAS_INLINE
1581     const_iterator2 find2 (int rank, size_type i, size_type j) const {
1582         const_subiterator2_type it2 (e_.find1 (rank, j, i));
1583 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1584         return const_iterator2 (*this, it2.index2 (), it2.index1 ());
1585 #else
1586         return const_iterator2 (*this, it2);
1587 #endif
1588     }
1589
1590     // Iterators enhance the iterators of the referenced expression
1591     // with the unary functor.
1592
1593 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1594     class const_iterator1:
1595         public container_const_reference<matrix_unary2>,
1596         public iterator_base_traits<typename E::const_iterator2::iterator_category>::template
1597         iterator_base<const_iterator1, value_type>::type {
1598     public:
1599         typedef typename E::const_iterator2::iterator_category iterator_category;
1600         typedef typename matrix_unary2::difference_type difference_type;
1601         typedef typename matrix_unary2::value_type value_type;
1602         typedef typename matrix_unary2::const_reference reference;
1603         typedef typename matrix_unary2::const_pointer pointer;
1604
1605         typedef const_iterator2 dual_iterator_type;
1606         typedef const_reverse_iterator2 dual_reverse_iterator_type;
1607
1608         // Construction and destruction
1609         BOOST_UBLAS_INLINE
1610         const_iterator1 ():
1611           container_const_reference<self_type> (), it_ () {}
1612         BOOST_UBLAS_INLINE
1613         const_iterator1 (const self_type &mu, const const_subiterator1_type &it):
1614           container_const_reference<self_type> (mu), it_ (it) {}
1615
1616         // Arithmetic
1617         BOOST_UBLAS_INLINE
1618         const_iterator1 &operator ++ () {
1619             ++ it_;
1620             return *this;
1621         }
1622         BOOST_UBLAS_INLINE
1623         const_iterator1 &operator -- () {
1624             -- it_;
1625             return *this;
1626         }
1627         BOOST_UBLAS_INLINE
1628         const_iterator1 &operator += (difference_type n) {
1629             it_ += n;
1630             return *this;
1631         }
1632         BOOST_UBLAS_INLINE
1633         const_iterator1 &operator -= (difference_type n) {
1634             it_ -= n;
1635             return *this;
1636         }
1637         BOOST_UBLAS_INLINE
1638         difference_type operator - (const const_iterator1 &it) const {
1639             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1640             return it_ - it.it_;
1641         }
1642
1643         // Dereference
1644         BOOST_UBLAS_INLINE
1645         const_reference operator * () const {
1646             return functor_type::apply (*it_);
1647         }
1648         BOOST_UBLAS_INLINE
1649         const_reference operator [] (difference_type n) const {
1650             return *(*this + n);
1651         }
1652
1653 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1654         BOOST_UBLAS_INLINE
1655 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1656         typename self_type::
1657 #endif
1658         const_iterator2 begin () const {
1659             return (*this) ().find2 (1, index1 (), 0);
1660         }
1661         BOOST_UBLAS_INLINE
1662 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1663         typename self_type::
1664 #endif
1665         const_iterator2 cbegin () const {
1666             return begin ();
1667         }
1668         BOOST_UBLAS_INLINE
1669 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1670         typename self_type::
1671 #endif
1672         const_iterator2 end () const {
1673             return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
1674         }
1675         BOOST_UBLAS_INLINE
1676 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1677         typename self_type::
1678 #endif
1679         const_iterator2 cend () const {
1680             return end ();
1681         }
1682         BOOST_UBLAS_INLINE
1683 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1684         typename self_type::
1685 #endif
1686         const_reverse_iterator2 rbegin () const {
1687             return const_reverse_iterator2 (end ());
1688         }
1689         BOOST_UBLAS_INLINE
1690 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1691         typename self_type::
1692 #endif
1693         const_reverse_iterator2 crbegin () const {
1694             return rbegin ();
1695         }
1696         BOOST_UBLAS_INLINE
1697 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1698         typename self_type::
1699 #endif
1700         const_reverse_iterator2 rend () const {
1701             return const_reverse_iterator2 (begin ());
1702         }
1703         BOOST_UBLAS_INLINE
1704 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1705         typename self_type::
1706 #endif
1707         const_reverse_iterator2 crend () const {
1708             return rend ();
1709         }
1710 #endif
1711
1712         // Indices
1713         BOOST_UBLAS_INLINE
1714         size_type index1 () const {
1715             return it_.index2 ();
1716         }
1717         BOOST_UBLAS_INLINE
1718         size_type index2 () const {
1719             return it_.index1 ();
1720         }
1721
1722         // Assignment
1723         BOOST_UBLAS_INLINE
1724         const_iterator1 &operator = (const const_iterator1 &it) {
1725             container_const_reference<self_type>::assign (&it ());
1726             it_ = it.it_;
1727             return *this;
1728         }
1729
1730         // Comparison
1731         BOOST_UBLAS_INLINE
1732         bool operator == (const const_iterator1 &it) const {
1733             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1734             return it_ == it.it_;
1735         }
1736         BOOST_UBLAS_INLINE
1737         bool operator < (const const_iterator1 &it) const {
1738             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1739             return it_ < it.it_;
1740         }
1741
1742     private:
1743         const_subiterator1_type it_;
1744     };
1745 #endif
1746
1747     BOOST_UBLAS_INLINE
1748     const_iterator1 begin1 () const {
1749         return find1 (0, 0, 0);
1750     }
1751     BOOST_UBLAS_INLINE
1752     const_iterator1 cbegin1 () const {
1753         return begin1 ();
1754     }
1755     BOOST_UBLAS_INLINE
1756     const_iterator1 end1 () const {
1757         return find1 (0, size1 (), 0);
1758     }
1759     BOOST_UBLAS_INLINE
1760     const_iterator1 cend1 () const {
1761         return end1 ();
1762     }
1763
1764 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1765     class const_iterator2:
1766         public container_const_reference<matrix_unary2>,
1767         public iterator_base_traits<typename E::const_iterator1::iterator_category>::template
1768         iterator_base<const_iterator2, value_type>::type {
1769     public:
1770         typedef typename E::const_iterator1::iterator_category iterator_category;
1771         typedef typename matrix_unary2::difference_type difference_type;
1772         typedef typename matrix_unary2::value_type value_type;
1773         typedef typename matrix_unary2::const_reference reference;
1774         typedef typename matrix_unary2::const_pointer pointer;
1775
1776         typedef const_iterator1 dual_iterator_type;
1777         typedef const_reverse_iterator1 dual_reverse_iterator_type;
1778
1779         // Construction and destruction
1780         BOOST_UBLAS_INLINE
1781         const_iterator2 ():
1782           container_const_reference<self_type> (), it_ () {}
1783         BOOST_UBLAS_INLINE
1784         const_iterator2 (const self_type &mu, const const_subiterator2_type &it):
1785           container_const_reference<self_type> (mu), it_ (it) {}
1786
1787         // Arithmetic
1788         BOOST_UBLAS_INLINE
1789         const_iterator2 &operator ++ () {
1790             ++ it_;
1791             return *this;
1792         }
1793         BOOST_UBLAS_INLINE
1794         const_iterator2 &operator -- () {
1795             -- it_;
1796             return *this;
1797         }
1798         BOOST_UBLAS_INLINE
1799         const_iterator2 &operator += (difference_type n) {
1800             it_ += n;
1801             return *this;
1802         }
1803         BOOST_UBLAS_INLINE
1804         const_iterator2 &operator -= (difference_type n) {
1805             it_ -= n;
1806             return *this;
1807         }
1808         BOOST_UBLAS_INLINE
1809         difference_type operator - (const const_iterator2 &it) const {
1810             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1811             return it_ - it.it_;
1812         }
1813
1814         // Dereference
1815         BOOST_UBLAS_INLINE
1816         const_reference operator * () const {
1817             return functor_type::apply (*it_);
1818         }
1819         BOOST_UBLAS_INLINE
1820         const_reference operator [] (difference_type n) const {
1821             return *(*this + n);
1822         }
1823
1824 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1825         BOOST_UBLAS_INLINE
1826 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1827         typename self_type::
1828 #endif
1829         const_iterator1 begin () const {
1830             return (*this) ().find1 (1, 0, index2 ());
1831         }
1832         BOOST_UBLAS_INLINE
1833 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1834         typename self_type::
1835 #endif
1836         const_iterator1 cbegin () const {
1837             return begin ();
1838         }
1839         BOOST_UBLAS_INLINE
1840 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1841         typename self_type::
1842 #endif
1843         const_iterator1 end () const {
1844             return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
1845         }
1846         BOOST_UBLAS_INLINE
1847 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1848         typename self_type::
1849 #endif
1850         const_iterator1 cend () const {
1851             return end ();
1852         }
1853         BOOST_UBLAS_INLINE
1854 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1855         typename self_type::
1856 #endif
1857         const_reverse_iterator1 rbegin () const {
1858             return const_reverse_iterator1 (end ());
1859         }
1860         BOOST_UBLAS_INLINE
1861 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1862         typename self_type::
1863 #endif
1864         const_reverse_iterator1 crbegin () const {
1865             return rbegin ();
1866         }
1867         BOOST_UBLAS_INLINE
1868 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1869         typename self_type::
1870 #endif
1871         const_reverse_iterator1 rend () const {
1872             return const_reverse_iterator1 (begin ());
1873         }
1874         BOOST_UBLAS_INLINE
1875 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1876         typename self_type::
1877 #endif
1878         const_reverse_iterator1 crend () const {
1879             return rend ();
1880         }
1881 #endif
1882
1883         // Indices
1884         BOOST_UBLAS_INLINE
1885         size_type index1 () const {
1886             return it_.index2 ();
1887         }
1888         BOOST_UBLAS_INLINE
1889         size_type index2 () const {
1890             return it_.index1 ();
1891         }
1892
1893         // Assignment
1894         BOOST_UBLAS_INLINE
1895         const_iterator2 &operator = (const const_iterator2 &it) {
1896             container_const_reference<self_type>::assign (&it ());
1897             it_ = it.it_;
1898             return *this;
1899         }
1900
1901         // Comparison
1902         BOOST_UBLAS_INLINE
1903         bool operator == (const const_iterator2 &it) const {
1904             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1905             return it_ == it.it_;
1906         }
1907         BOOST_UBLAS_INLINE
1908         bool operator < (const const_iterator2 &it) const {
1909             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1910             return it_ < it.it_;
1911         }
1912
1913     private:
1914         const_subiterator2_type it_;
1915     };
1916 #endif
1917
1918     BOOST_UBLAS_INLINE
1919     const_iterator2 begin2 () const {
1920         return find2 (0, 0, 0);
1921     }
1922     BOOST_UBLAS_INLINE
1923     const_iterator2 cbegin2 () const {
1924         return begin2 ();
1925     }
1926     BOOST_UBLAS_INLINE
1927     const_iterator2 end2 () const {
1928         return find2 (0, 0, size2 ());
1929     }
1930     BOOST_UBLAS_INLINE
1931     const_iterator2 cend2 () const {
1932         return end2 ();
1933     }
1934
1935     // Reverse iterators
1936
1937     BOOST_UBLAS_INLINE
1938     const_reverse_iterator1 rbegin1 () const {
1939         return const_reverse_iterator1 (end1 ());
1940     }
1941     BOOST_UBLAS_INLINE
1942     const_reverse_iterator1 crbegin1 () const {
1943         return rbegin1 ();
1944     }
1945     BOOST_UBLAS_INLINE
1946     const_reverse_iterator1 rend1 () const {
1947         return const_reverse_iterator1 (begin1 ());
1948     }
1949     BOOST_UBLAS_INLINE
1950     const_reverse_iterator1 crend1 () const {
1951         return rend1 ();
1952     }
1953
1954     BOOST_UBLAS_INLINE
1955     const_reverse_iterator2 rbegin2 () const {
1956         return const_reverse_iterator2 (end2 ());
1957     }
1958     BOOST_UBLAS_INLINE
1959     const_reverse_iterator2 crbegin2 () const {
1960         return rbegin2 ();
1961     }
1962     BOOST_UBLAS_INLINE
1963     const_reverse_iterator2 rend2 () const {
1964         return const_reverse_iterator2 (begin2 ());
1965     }
1966     BOOST_UBLAS_INLINE
1967     const_reverse_iterator2 crend2 () const {
1968         return rend2 ();
1969     }
1970
1971 private:
1972     expression_closure_type e_;
1973 };
1974
1975 template<class E, class F>
1976 struct matrix_unary2_traits {
1977     typedef matrix_unary2<E, F> expression_type;
1978 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
1979     typedef expression_type result_type;
1980 #else
1981     typedef typename E::matrix_temporary_type result_type;
1982 #endif
1983 };
1984
1985 // (trans m) [i] [j] = m [j] [i]
1986 template<class E>
1987 BOOST_UBLAS_INLINE
1988 typename matrix_unary2_traits<const E, scalar_identity<typename E::value_type> >::result_type
1989 trans (const matrix_expression<E> &e) {
1990     typedef typename matrix_unary2_traits<const E, scalar_identity<typename E::value_type> >::expression_type expression_type;
1991     return expression_type (e ());
1992 }
1993 template<class E>
1994 BOOST_UBLAS_INLINE
1995 typename matrix_unary2_traits<E, scalar_identity<typename E::value_type> >::result_type
1996 trans (matrix_expression<E> &e) {
1997     typedef typename matrix_unary2_traits<E, scalar_identity<typename E::value_type> >::expression_type expression_type;
1998     return expression_type (e ());
1999 }
2000
2001 // (herm m) [i] [j] = conj (m [j] [i])
2002 template<class E>
2003 BOOST_UBLAS_INLINE
2004 typename matrix_unary2_traits<E, scalar_conj<typename E::value_type> >::result_type
2005 herm (const matrix_expression<E> &e) {
2006     typedef typename matrix_unary2_traits<E, scalar_conj<typename E::value_type> >::expression_type expression_type;
2007     return expression_type (e ());
2008 }
2009
2010 template<class E1, class E2, class F>
2011 class matrix_binary:
2012     public matrix_expression<matrix_binary<E1, E2, F> > {
2013
2014     typedef E1 expression1_type;
2015     typedef E2 expression2_type;
2016     typedef F functor_type;
2017 public:
2018     typedef typename E1::const_closure_type expression1_closure_type;
2019     typedef typename E2::const_closure_type expression2_closure_type;
2020 private:
2021     typedef matrix_binary<E1, E2, F> self_type;
2022 public:
2023 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
2024     using matrix_expression<self_type>::operator ();
2025 #endif
2026     typedef typename promote_traits<typename E1::size_type, typename E2::size_type>::promote_type size_type;
2027     typedef typename promote_traits<typename E1::difference_type, typename E2::difference_type>::promote_type difference_type;
2028     typedef typename F::result_type value_type;
2029     typedef value_type const_reference;
2030     typedef const_reference reference;
2031     typedef const self_type const_closure_type;
2032     typedef const_closure_type closure_type;
2033     typedef unknown_orientation_tag orientation_category;
2034     typedef unknown_storage_tag storage_category;
2035
2036     // Construction and destruction
2037     BOOST_UBLAS_INLINE
2038     matrix_binary (const E1 &e1, const E2 &e2):
2039       e1_ (e1), e2_ (e2) {}
2040
2041     // Accessors
2042     BOOST_UBLAS_INLINE
2043     size_type size1 () const {
2044         return BOOST_UBLAS_SAME (e1_.size1 (), e2_.size1 ());
2045     }
2046     BOOST_UBLAS_INLINE
2047     size_type size2 () const {
2048         return BOOST_UBLAS_SAME (e1_.size2 (), e2_.size2 ());
2049     }
2050
2051 public:
2052     // Expression accessors
2053     BOOST_UBLAS_INLINE
2054     const expression1_closure_type &expression1 () const {
2055         return e1_;
2056     }
2057     BOOST_UBLAS_INLINE
2058     const expression2_closure_type &expression2 () const {
2059         return e2_;
2060     }
2061
2062 public:
2063     // Element access
2064     BOOST_UBLAS_INLINE
2065     const_reference operator () (size_type i, size_type j) const {
2066         return functor_type::apply (e1_ (i, j), e2_ (i, j));
2067     }
2068
2069     // Element access
2070     BOOST_UBLAS_INLINE
2071     const_reference operator () (size_type i) const {
2072         return functor_type::apply (e1_ (i), e2_ (i));
2073     }
2074
2075     // Closure comparison
2076     BOOST_UBLAS_INLINE
2077     bool same_closure (const matrix_binary &mb) const {
2078         return (*this).expression1 ().same_closure (mb.expression1 ()) &&
2079             (*this).expression2 ().same_closure (mb.expression2 ());
2080     }
2081
2082     // Iterator types
2083 private:
2084     typedef typename E1::const_iterator1 const_iterator11_type;
2085     typedef typename E1::const_iterator2 const_iterator12_type;
2086     typedef typename E2::const_iterator1 const_iterator21_type;
2087     typedef typename E2::const_iterator2 const_iterator22_type;
2088     typedef const value_type *const_pointer;
2089
2090 public:
2091 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
2092     typedef typename iterator_restrict_traits<typename const_iterator11_type::iterator_category,
2093     typename const_iterator21_type::iterator_category>::iterator_category iterator_category1;
2094     typedef indexed_const_iterator1<const_closure_type, iterator_category1> const_iterator1;
2095     typedef const_iterator1 iterator1;
2096     typedef typename iterator_restrict_traits<typename const_iterator12_type::iterator_category,
2097     typename const_iterator22_type::iterator_category>::iterator_category iterator_category2;
2098     typedef indexed_const_iterator2<const_closure_type, iterator_category2> const_iterator2;
2099     typedef const_iterator2 iterator2;
2100 #else
2101     class const_iterator1;
2102     typedef const_iterator1 iterator1;
2103     class const_iterator2;
2104     typedef const_iterator2 iterator2;
2105 #endif
2106     typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
2107     typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
2108
2109     // Element lookup
2110     BOOST_UBLAS_INLINE
2111     const_iterator1 find1 (int rank, size_type i, size_type j) const {
2112         const_iterator11_type it11 (e1_.find1 (rank, i, j));
2113         const_iterator11_type it11_end (e1_.find1 (rank, size1 (), j));
2114         const_iterator21_type it21 (e2_.find1 (rank, i, j));
2115         const_iterator21_type it21_end (e2_.find1 (rank, size1 (), j));
2116         BOOST_UBLAS_CHECK (rank == 0 || it11 == it11_end || it11.index2 () == j, internal_logic ())
2117             BOOST_UBLAS_CHECK (rank == 0 || it21 == it21_end || it21.index2 () == j, internal_logic ())
2118             i = (std::min) (it11 != it11_end ? it11.index1 () : size1 (),
2119                             it21 != it21_end ? it21.index1 () : size1 ());
2120 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
2121         return const_iterator1 (*this, i, j);
2122 #else
2123         return const_iterator1 (*this, i, j, it11, it11_end, it21, it21_end);
2124 #endif
2125     }
2126     BOOST_UBLAS_INLINE
2127     const_iterator2 find2 (int rank, size_type i, size_type j) const {
2128         const_iterator12_type it12 (e1_.find2 (rank, i, j));
2129         const_iterator12_type it12_end (e1_.find2 (rank, i, size2 ()));
2130         const_iterator22_type it22 (e2_.find2 (rank, i, j));
2131         const_iterator22_type it22_end (e2_.find2 (rank, i, size2 ()));
2132         BOOST_UBLAS_CHECK (rank == 0 || it12 == it12_end || it12.index1 () == i, internal_logic ())
2133             BOOST_UBLAS_CHECK (rank == 0 || it22 == it22_end || it22.index1 () == i, internal_logic ())
2134             j = (std::min) (it12 != it12_end ? it12.index2 () : size2 (),
2135                             it22 != it22_end ? it22.index2 () : size2 ());
2136 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
2137         return const_iterator2 (*this, i, j);
2138 #else
2139         return const_iterator2 (*this, i, j, it12, it12_end, it22, it22_end);
2140 #endif
2141     }
2142
2143     // Iterators enhance the iterators of the referenced expression
2144     // with the binary functor.
2145
2146 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
2147     class const_iterator1:
2148         public container_const_reference<matrix_binary>,
2149         public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
2150         typename E2::const_iterator1::iterator_category>::iterator_category>::template
2151         iterator_base<const_iterator1, value_type>::type {
2152     public:
2153         typedef typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
2154         typename E2::const_iterator1::iterator_category>::iterator_category iterator_category;
2155         typedef typename matrix_binary::difference_type difference_type;
2156         typedef typename matrix_binary::value_type value_type;
2157         typedef typename matrix_binary::const_reference reference;
2158         typedef typename matrix_binary::const_pointer pointer;
2159
2160         typedef const_iterator2 dual_iterator_type;
2161         typedef const_reverse_iterator2 dual_reverse_iterator_type;
2162
2163         // Construction and destruction
2164         BOOST_UBLAS_INLINE
2165         const_iterator1 ():
2166           container_const_reference<self_type> (), i_ (), j_ (), it1_ (), it1_end_ (), it2_ (), it2_end_ () {}
2167         BOOST_UBLAS_INLINE
2168         const_iterator1 (const self_type &mb, size_type i, size_type j,
2169                          const const_iterator11_type &it1, const const_iterator11_type &it1_end,
2170                          const const_iterator21_type &it2, const const_iterator21_type &it2_end):
2171           container_const_reference<self_type> (mb), i_ (i), j_ (j), it1_ (it1), it1_end_ (it1_end), it2_ (it2), it2_end_ (it2_end) {}
2172
2173     private:
2174         // Dense specializations
2175         BOOST_UBLAS_INLINE
2176         void increment (dense_random_access_iterator_tag) {
2177             ++ i_; ++ it1_; ++ it2_;
2178         }
2179         BOOST_UBLAS_INLINE
2180         void decrement (dense_random_access_iterator_tag) {
2181             -- i_; -- it1_; -- it2_;
2182         }
2183         BOOST_UBLAS_INLINE
2184         void increment (dense_random_access_iterator_tag, difference_type n) {
2185             i_ += n; it1_ += n; it2_ += n;
2186         }
2187         BOOST_UBLAS_INLINE
2188         void decrement (dense_random_access_iterator_tag, difference_type n) {
2189             i_ -= n; it1_ -= n; it2_ -= n;
2190         }
2191         BOOST_UBLAS_INLINE
2192         value_type dereference (dense_random_access_iterator_tag) const {
2193             return functor_type::apply (*it1_, *it2_);
2194         }
2195
2196         // Packed specializations
2197         BOOST_UBLAS_INLINE
2198         void increment (packed_random_access_iterator_tag) {
2199             if (it1_ != it1_end_)
2200                 if (it1_.index1 () <= i_)
2201                     ++ it1_;
2202             if (it2_ != it2_end_)
2203                 if (it2_.index1 () <= i_)
2204                     ++ it2_;
2205             ++ i_;
2206         }
2207         BOOST_UBLAS_INLINE
2208         void decrement (packed_random_access_iterator_tag) {
2209             if (it1_ != it1_end_)
2210                 if (i_ <= it1_.index1 ())
2211                     -- it1_;
2212             if (it2_ != it2_end_)
2213                 if (i_ <= it2_.index1 ())
2214                     -- it2_;
2215             -- i_;
2216         }
2217         BOOST_UBLAS_INLINE
2218         void increment (packed_random_access_iterator_tag, difference_type n) {
2219             while (n > 0) {
2220                 increment (packed_random_access_iterator_tag ());
2221                 --n;
2222             }
2223             while (n < 0) {
2224                 decrement (packed_random_access_iterator_tag ());
2225                 ++n;
2226             }
2227         }
2228         BOOST_UBLAS_INLINE
2229         void decrement (packed_random_access_iterator_tag, difference_type n) {
2230             while (n > 0) {
2231                 decrement (packed_random_access_iterator_tag ());
2232                 --n;
2233             }
2234             while (n < 0) {
2235                 increment (packed_random_access_iterator_tag ());
2236                 ++n;
2237             }
2238         }
2239         BOOST_UBLAS_INLINE
2240         value_type dereference (packed_random_access_iterator_tag) const {
2241             typename E1::value_type t1 = typename E1::value_type/*zero*/();
2242             if (it1_ != it1_end_) {
2243                 BOOST_UBLAS_CHECK (it1_.index2 () == j_, internal_logic ());
2244                 if (it1_.index1 () == i_)
2245                     t1 = *it1_;
2246             }
2247             typename E2::value_type t2 = typename E2::value_type/*zero*/();
2248             if (it2_ != it2_end_) {
2249                 BOOST_UBLAS_CHECK (it2_.index2 () == j_, internal_logic ());
2250                 if (it2_.index1 () == i_)
2251                     t2 = *it2_;
2252             }
2253             return functor_type::apply (t1, t2);
2254         }
2255
2256         // Sparse specializations
2257         BOOST_UBLAS_INLINE
2258         void increment (sparse_bidirectional_iterator_tag) {
2259             size_type index1 = (*this) ().size1 ();
2260             if (it1_ != it1_end_) {
2261                 if (it1_.index1 () <= i_)
2262                     ++ it1_;
2263                 if (it1_ != it1_end_)
2264                     index1 = it1_.index1 ();
2265             }
2266             size_type index2 = (*this) ().size1 ();
2267             if (it2_ != it2_end_)
2268                 if (it2_.index1 () <= i_)
2269                     ++ it2_;
2270             if (it2_ != it2_end_) {
2271                 index2 = it2_.index1 ();
2272             }
2273             i_ = (std::min) (index1, index2);
2274         }
2275         BOOST_UBLAS_INLINE
2276         void decrement (sparse_bidirectional_iterator_tag) {
2277             size_type index1 = (*this) ().size1 ();
2278             if (it1_ != it1_end_) {
2279                 if (i_ <= it1_.index1 ())
2280                     -- it1_;
2281                 if (it1_ != it1_end_)
2282                     index1 = it1_.index1 ();
2283             }
2284             size_type index2 = (*this) ().size1 ();
2285             if (it2_ != it2_end_) {
2286                 if (i_ <= it2_.index1 ())
2287                     -- it2_;
2288                 if (it2_ != it2_end_)
2289                     index2 = it2_.index1 ();
2290             }
2291             i_ = (std::max) (index1, index2);
2292         }
2293         BOOST_UBLAS_INLINE
2294         void increment (sparse_bidirectional_iterator_tag, difference_type n) {
2295             while (n > 0) {
2296                 increment (sparse_bidirectional_iterator_tag ());
2297                 --n;
2298             }
2299             while (n < 0) {
2300                 decrement (sparse_bidirectional_iterator_tag ());
2301                 ++n;
2302             }
2303         }
2304         BOOST_UBLAS_INLINE
2305         void decrement (sparse_bidirectional_iterator_tag, difference_type n) {
2306             while (n > 0) {
2307                 decrement (sparse_bidirectional_iterator_tag ());
2308                 --n;
2309             }
2310             while (n < 0) {
2311                 increment (sparse_bidirectional_iterator_tag ());
2312                 ++n;
2313             }
2314         }
2315         BOOST_UBLAS_INLINE
2316         value_type dereference (sparse_bidirectional_iterator_tag) const {
2317             typename E1::value_type t1 = typename E1::value_type/*zero*/();
2318             if (it1_ != it1_end_) {
2319                 BOOST_UBLAS_CHECK (it1_.index2 () == j_, internal_logic ());
2320                 if (it1_.index1 () == i_)
2321                     t1 = *it1_;
2322             }
2323             typename E2::value_type t2 = typename E2::value_type/*zero*/();
2324             if (it2_ != it2_end_) {
2325                 BOOST_UBLAS_CHECK (it2_.index2 () == j_, internal_logic ());
2326                 if (it2_.index1 () == i_)
2327                     t2 = *it2_;
2328             }
2329             return functor_type::apply (t1, t2);
2330         }
2331
2332     public:
2333         // Arithmetic
2334         BOOST_UBLAS_INLINE
2335         const_iterator1 &operator ++ () {
2336             increment (iterator_category ());
2337             return *this;
2338         }
2339         BOOST_UBLAS_INLINE
2340         const_iterator1 &operator -- () {
2341             decrement (iterator_category ());
2342             return *this;
2343         }
2344         BOOST_UBLAS_INLINE
2345         const_iterator1 &operator += (difference_type n) {
2346             increment (iterator_category (), n);
2347             return *this;
2348         }
2349         BOOST_UBLAS_INLINE
2350         const_iterator1 &operator -= (difference_type n) {
2351             decrement (iterator_category (), n);
2352             return *this;
2353         }
2354         BOOST_UBLAS_INLINE
2355         difference_type operator - (const const_iterator1 &it) const {
2356             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2357             BOOST_UBLAS_CHECK (index2 () == it.index2 (), external_logic ());
2358             return index1 () - it.index1 ();
2359         }
2360
2361         // Dereference
2362         BOOST_UBLAS_INLINE
2363         const_reference operator * () const {
2364             return dereference (iterator_category ());
2365         }
2366         BOOST_UBLAS_INLINE
2367         const_reference operator [] (difference_type n) const {
2368             return *(*this + n);
2369         }
2370
2371 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
2372         BOOST_UBLAS_INLINE
2373 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2374         typename self_type::
2375 #endif
2376         const_iterator2 begin () const {
2377             return (*this) ().find2 (1, index1 (), 0);
2378         }
2379         BOOST_UBLAS_INLINE
2380 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2381         typename self_type::
2382 #endif
2383         const_iterator2 cbegin () const {
2384             return begin ();
2385         }
2386         BOOST_UBLAS_INLINE
2387 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2388         typename self_type::
2389 #endif
2390         const_iterator2 end () const {
2391             return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
2392         }
2393         BOOST_UBLAS_INLINE
2394 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2395         typename self_type::
2396 #endif
2397         const_iterator2 cend () const {
2398             return end ();
2399         }
2400         BOOST_UBLAS_INLINE
2401 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2402         typename self_type::
2403 #endif
2404         const_reverse_iterator2 rbegin () const {
2405             return const_reverse_iterator2 (end ());
2406         }
2407         BOOST_UBLAS_INLINE
2408 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2409         typename self_type::
2410 #endif
2411         const_reverse_iterator2 crbegin () const {
2412             return rbegin ();
2413         }
2414         BOOST_UBLAS_INLINE
2415 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2416         typename self_type::
2417 #endif
2418         const_reverse_iterator2 rend () const {
2419             return const_reverse_iterator2 (begin ());
2420         }
2421         BOOST_UBLAS_INLINE
2422 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2423         typename self_type::
2424 #endif
2425         const_reverse_iterator2 crend () const {
2426             return rend ();
2427         }
2428 #endif
2429
2430         // Indices
2431         BOOST_UBLAS_INLINE
2432         size_type index1 () const {
2433             return i_;
2434         }
2435         BOOST_UBLAS_INLINE
2436         size_type index2 () const {
2437             // if (it1_ != it1_end_ && it2_ != it2_end_)
2438             //    return BOOST_UBLAS_SAME (it1_.index2 (), it2_.index2 ());
2439             // else
2440             return j_;
2441         }
2442
2443         // Assignment
2444         BOOST_UBLAS_INLINE
2445         const_iterator1 &operator = (const const_iterator1 &it) {
2446             container_const_reference<self_type>::assign (&it ());
2447             i_ = it.i_;
2448             j_ = it.j_;
2449             it1_ = it.it1_;
2450             it1_end_ = it.it1_end_;
2451             it2_ = it.it2_;
2452             it2_end_ = it.it2_end_;
2453             return *this;
2454         }
2455
2456         // Comparison
2457         BOOST_UBLAS_INLINE
2458         bool operator == (const const_iterator1 &it) const {
2459             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2460             BOOST_UBLAS_CHECK (index2 () == it.index2 (), external_logic ());
2461             return index1 () == it.index1 ();
2462         }
2463         BOOST_UBLAS_INLINE
2464         bool operator < (const const_iterator1 &it) const {
2465             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2466             BOOST_UBLAS_CHECK (index2 () == it.index2 (), external_logic ());
2467             return index1 () < it.index1 ();
2468         }
2469
2470     private:
2471         size_type i_;
2472         size_type j_;
2473         const_iterator11_type it1_;
2474         const_iterator11_type it1_end_;
2475         const_iterator21_type it2_;
2476         const_iterator21_type it2_end_;
2477     };
2478 #endif
2479
2480     BOOST_UBLAS_INLINE
2481     const_iterator1 begin1 () const {
2482         return find1 (0, 0, 0);
2483     }
2484     BOOST_UBLAS_INLINE
2485     const_iterator1 cbegin1 () const {
2486         return begin1 ();
2487     }
2488     BOOST_UBLAS_INLINE
2489     const_iterator1 end1 () const {
2490         return find1 (0, size1 (), 0);
2491     }
2492     BOOST_UBLAS_INLINE
2493     const_iterator1 cend1 () const {
2494         return end1 ();
2495     }
2496
2497 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
2498     class const_iterator2:
2499         public container_const_reference<matrix_binary>,
2500         public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator2::iterator_category,
2501         typename E2::const_iterator2::iterator_category>::iterator_category>::template
2502         iterator_base<const_iterator2, value_type>::type {
2503     public:
2504         typedef typename iterator_restrict_traits<typename E1::const_iterator2::iterator_category,
2505         typename E2::const_iterator2::iterator_category>::iterator_category iterator_category;
2506         typedef typename matrix_binary::difference_type difference_type;
2507         typedef typename matrix_binary::value_type value_type;
2508         typedef typename matrix_binary::const_reference reference;
2509         typedef typename matrix_binary::const_pointer pointer;
2510
2511         typedef const_iterator1 dual_iterator_type;
2512         typedef const_reverse_iterator1 dual_reverse_iterator_type;
2513
2514         // Construction and destruction
2515         BOOST_UBLAS_INLINE
2516         const_iterator2 ():
2517           container_const_reference<self_type> (), i_ (), j_ (), it1_ (), it1_end_ (), it2_ (), it2_end_ () {}
2518         BOOST_UBLAS_INLINE
2519         const_iterator2 (const self_type &mb, size_type i, size_type j,
2520                          const const_iterator12_type &it1, const const_iterator12_type &it1_end,
2521                          const const_iterator22_type &it2, const const_iterator22_type &it2_end):
2522           container_const_reference<self_type> (mb), i_ (i), j_ (j), it1_ (it1), it1_end_ (it1_end), it2_ (it2), it2_end_ (it2_end) {}
2523
2524     private:
2525         // Dense access specializations
2526         BOOST_UBLAS_INLINE
2527         void increment (dense_random_access_iterator_tag) {
2528             ++ j_; ++ it1_; ++ it2_;
2529         }
2530         BOOST_UBLAS_INLINE
2531         void decrement (dense_random_access_iterator_tag) {
2532             -- j_; -- it1_; -- it2_;
2533         }
2534         BOOST_UBLAS_INLINE
2535         void increment (dense_random_access_iterator_tag, difference_type n) {
2536             j_ += n; it1_ += n; it2_ += n;
2537         }
2538         BOOST_UBLAS_INLINE
2539         void decrement (dense_random_access_iterator_tag, difference_type n) {
2540             j_ -= n; it1_ -= n; it2_ -= n;
2541         }
2542         BOOST_UBLAS_INLINE
2543         value_type dereference (dense_random_access_iterator_tag) const {
2544             return functor_type::apply (*it1_, *it2_);
2545         }
2546
2547         // Packed specializations
2548         BOOST_UBLAS_INLINE
2549         void increment (packed_random_access_iterator_tag) {
2550             if (it1_ != it1_end_)
2551                 if (it1_.index2 () <= j_)
2552                     ++ it1_;
2553             if (it2_ != it2_end_)
2554                 if (it2_.index2 () <= j_)
2555                     ++ it2_;
2556             ++ j_;
2557         }
2558         BOOST_UBLAS_INLINE
2559         void decrement (packed_random_access_iterator_tag) {
2560             if (it1_ != it1_end_)
2561                 if (j_ <= it1_.index2 ())
2562                     -- it1_;
2563             if (it2_ != it2_end_)
2564                 if (j_ <= it2_.index2 ())
2565                     -- it2_;
2566             -- j_;
2567         }
2568         BOOST_UBLAS_INLINE
2569         void increment (packed_random_access_iterator_tag, difference_type n) {
2570             while (n > 0) {
2571                 increment (packed_random_access_iterator_tag ());
2572                 --n;
2573             }
2574             while (n < 0) {
2575                 decrement (packed_random_access_iterator_tag ());
2576                 ++n;
2577             }
2578         }
2579         BOOST_UBLAS_INLINE
2580         void decrement (packed_random_access_iterator_tag, difference_type n) {
2581             while (n > 0) {
2582                 decrement (packed_random_access_iterator_tag ());
2583                 --n;
2584             }
2585             while (n < 0) {
2586                 increment (packed_random_access_iterator_tag ());
2587                 ++n;
2588             }
2589         }
2590         BOOST_UBLAS_INLINE
2591         value_type dereference (packed_random_access_iterator_tag) const {
2592             typename E1::value_type t1 = typename E1::value_type/*zero*/();
2593             if (it1_ != it1_end_) {
2594                 BOOST_UBLAS_CHECK (it1_.index1 () == i_, internal_logic ());
2595                 if (it1_.index2 () == j_)
2596                     t1 = *it1_;
2597             }
2598             typename E2::value_type t2 = typename E2::value_type/*zero*/();
2599             if (it2_ != it2_end_) {
2600                 BOOST_UBLAS_CHECK (it2_.index1 () == i_, internal_logic ());
2601                 if (it2_.index2 () == j_)
2602                     t2 = *it2_;
2603             }
2604             return functor_type::apply (t1, t2);
2605         }
2606
2607         // Sparse specializations
2608         BOOST_UBLAS_INLINE
2609         void increment (sparse_bidirectional_iterator_tag) {
2610             size_type index1 = (*this) ().size2 ();
2611             if (it1_ != it1_end_) {
2612                 if (it1_.index2 () <= j_)
2613                     ++ it1_;
2614                 if (it1_ != it1_end_)
2615                     index1 = it1_.index2 ();
2616             }
2617             size_type index2 = (*this) ().size2 ();
2618             if (it2_ != it2_end_) {
2619                 if (it2_.index2 () <= j_)
2620                     ++ it2_;
2621                 if (it2_ != it2_end_)
2622                     index2 = it2_.index2 ();
2623             }
2624             j_ = (std::min) (index1, index2);
2625         }
2626         BOOST_UBLAS_INLINE
2627         void decrement (sparse_bidirectional_iterator_tag) {
2628             size_type index1 = (*this) ().size2 ();
2629             if (it1_ != it1_end_) {
2630                 if (j_ <= it1_.index2 ())
2631                     -- it1_;
2632                 if (it1_ != it1_end_)
2633                     index1 = it1_.index2 ();
2634             }
2635             size_type index2 = (*this) ().size2 ();
2636             if (it2_ != it2_end_) {
2637                 if (j_ <= it2_.index2 ())
2638                     -- it2_;
2639                 if (it2_ != it2_end_)
2640                     index2 = it2_.index2 ();
2641             }
2642             j_ = (std::max) (index1, index2);
2643         }
2644         BOOST_UBLAS_INLINE
2645         void increment (sparse_bidirectional_iterator_tag, difference_type n) {
2646             while (n > 0) {
2647                 increment (sparse_bidirectional_iterator_tag ());
2648                 --n;
2649             }
2650             while (n < 0) {
2651                 decrement (sparse_bidirectional_iterator_tag ());
2652                 ++n;
2653             }
2654         }
2655         BOOST_UBLAS_INLINE
2656         void decrement (sparse_bidirectional_iterator_tag, difference_type n) {
2657             while (n > 0) {
2658                 decrement (sparse_bidirectional_iterator_tag ());
2659                 --n;
2660             }
2661             while (n < 0) {
2662                 increment (sparse_bidirectional_iterator_tag ());
2663                 ++n;
2664             }
2665         }
2666         BOOST_UBLAS_INLINE
2667         value_type dereference (sparse_bidirectional_iterator_tag) const {
2668             typename E1::value_type t1 = typename E1::value_type/*zero*/();
2669             if (it1_ != it1_end_) {
2670                 BOOST_UBLAS_CHECK (it1_.index1 () == i_, internal_logic ());
2671                 if (it1_.index2 () == j_)
2672                     t1 = *it1_;
2673             }
2674             typename E2::value_type t2 = typename E2::value_type/*zero*/();
2675             if (it2_ != it2_end_) {
2676                 BOOST_UBLAS_CHECK (it2_.index1 () == i_, internal_logic ());
2677                 if (it2_.index2 () == j_)
2678                     t2 = *it2_;
2679             }
2680             return functor_type::apply (t1, t2);
2681         }
2682
2683     public:
2684         // Arithmetic
2685         BOOST_UBLAS_INLINE
2686         const_iterator2 &operator ++ () {
2687             increment (iterator_category ());
2688             return *this;
2689         }
2690         BOOST_UBLAS_INLINE
2691         const_iterator2 &operator -- () {
2692             decrement (iterator_category ());
2693             return *this;
2694         }
2695         BOOST_UBLAS_INLINE
2696         const_iterator2 &operator += (difference_type n) {
2697             increment (iterator_category (), n);
2698             return *this;
2699         }
2700         BOOST_UBLAS_INLINE
2701         const_iterator2 &operator -= (difference_type n) {
2702             decrement (iterator_category (), n);
2703             return *this;
2704         }
2705         BOOST_UBLAS_INLINE
2706         difference_type operator - (const const_iterator2 &it) const {
2707             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2708             BOOST_UBLAS_CHECK (index1 () == it.index1 (), external_logic ());
2709             return index2 () - it.index2 ();
2710         }
2711
2712         // Dereference
2713         BOOST_UBLAS_INLINE
2714         const_reference operator * () const {
2715             return dereference (iterator_category ());
2716         }
2717         BOOST_UBLAS_INLINE
2718         const_reference operator [] (difference_type n) const {
2719             return *(*this + n);
2720         }
2721
2722 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
2723         BOOST_UBLAS_INLINE
2724 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2725         typename self_type::
2726 #endif
2727         const_iterator1 begin () const {
2728             return (*this) ().find1 (1, 0, index2 ());
2729         }
2730         BOOST_UBLAS_INLINE
2731 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2732         typename self_type::
2733 #endif
2734         const_iterator1 cbegin () const {
2735             return begin ();
2736         }
2737         BOOST_UBLAS_INLINE
2738 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2739         typename self_type::
2740 #endif
2741         const_iterator1 end () const {
2742             return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
2743         }
2744         BOOST_UBLAS_INLINE
2745 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2746         typename self_type::
2747 #endif
2748         const_iterator1 cend () const {
2749             return end ();
2750         }
2751         BOOST_UBLAS_INLINE
2752 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2753         typename self_type::
2754 #endif
2755         const_reverse_iterator1 rbegin () const {
2756             return const_reverse_iterator1 (end ());
2757         }
2758         BOOST_UBLAS_INLINE
2759 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2760         typename self_type::
2761 #endif
2762         const_reverse_iterator1 crbegin () const {
2763             return rbegin ();
2764         }
2765         BOOST_UBLAS_INLINE
2766 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2767         typename self_type::
2768 #endif
2769         const_reverse_iterator1 rend () const {
2770             return const_reverse_iterator1 (begin ());
2771         }
2772         BOOST_UBLAS_INLINE
2773 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2774         typename self_type::
2775 #endif
2776         const_reverse_iterator1 crend () const {
2777             return rend ();
2778         }
2779 #endif
2780
2781         // Indices
2782         BOOST_UBLAS_INLINE
2783         size_type index1 () const {
2784             // if (it1_ != it1_end_ && it2_ != it2_end_)
2785             //    return BOOST_UBLAS_SAME (it1_.index1 (), it2_.index1 ());
2786             // else
2787             return i_;
2788         }
2789         BOOST_UBLAS_INLINE
2790         size_type index2 () const {
2791             return j_;
2792         }
2793
2794         // Assignment
2795         BOOST_UBLAS_INLINE
2796         const_iterator2 &operator = (const const_iterator2 &it) {
2797             container_const_reference<self_type>::assign (&it ());
2798             i_ = it.i_;
2799             j_ = it.j_;
2800             it1_ = it.it1_;
2801             it1_end_ = it.it1_end_;
2802             it2_ = it.it2_;
2803             it2_end_ = it.it2_end_;
2804             return *this;
2805         }
2806
2807         // Comparison
2808         BOOST_UBLAS_INLINE
2809         bool operator == (const const_iterator2 &it) const {
2810             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2811             BOOST_UBLAS_CHECK (index1 () == it.index1 (), external_logic ());
2812             return index2 () == it.index2 ();
2813         }
2814         BOOST_UBLAS_INLINE
2815         bool operator < (const const_iterator2 &it) const {
2816             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2817             BOOST_UBLAS_CHECK (index1 () == it.index1 (), external_logic ());
2818             return index2 () < it.index2 ();
2819         }
2820
2821     private:
2822         size_type i_;
2823         size_type j_;
2824         const_iterator12_type it1_;
2825         const_iterator12_type it1_end_;
2826         const_iterator22_type it2_;
2827         const_iterator22_type it2_end_;
2828     };
2829 #endif
2830
2831     BOOST_UBLAS_INLINE
2832     const_iterator2 begin2 () const {
2833         return find2 (0, 0, 0);
2834     }
2835     BOOST_UBLAS_INLINE
2836     const_iterator2 cbegin2 () const {
2837         return begin2 ();
2838     }
2839     BOOST_UBLAS_INLINE
2840     const_iterator2 end2 () const {
2841         return find2 (0, 0, size2 ());
2842     }
2843     BOOST_UBLAS_INLINE
2844     const_iterator2 cend2 () const {
2845         return end2 ();
2846     }
2847
2848     // Reverse iterators
2849
2850     BOOST_UBLAS_INLINE
2851     const_reverse_iterator1 rbegin1 () const {
2852         return const_reverse_iterator1 (end1 ());
2853     }
2854     BOOST_UBLAS_INLINE
2855     const_reverse_iterator1 crbegin1 () const {
2856         return rbegin1 ();
2857     }
2858     BOOST_UBLAS_INLINE
2859     const_reverse_iterator1 rend1 () const {
2860         return const_reverse_iterator1 (begin1 ());
2861     }
2862     BOOST_UBLAS_INLINE
2863     const_reverse_iterator1 crend1 () const {
2864         return rend1 ();
2865     }
2866
2867     BOOST_UBLAS_INLINE
2868     const_reverse_iterator2 rbegin2 () const {
2869         return const_reverse_iterator2 (end2 ());
2870     }
2871     BOOST_UBLAS_INLINE
2872     const_reverse_iterator2 crbegin2 () const {
2873         return rbegin2 ();
2874     }
2875     BOOST_UBLAS_INLINE
2876     const_reverse_iterator2 rend2 () const {
2877         return const_reverse_iterator2 (begin2 ());
2878     }
2879     BOOST_UBLAS_INLINE
2880     const_reverse_iterator2 crend2 () const {
2881         return rend2 ();
2882     }
2883
2884 private:
2885     expression1_closure_type e1_;
2886     expression2_closure_type e2_;
2887 };
2888
2889 template<class E1, class E2, class F>
2890 struct matrix_binary_traits {
2891     typedef matrix_binary<E1, E2, F> expression_type;
2892 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
2893     typedef expression_type result_type;
2894 #else
2895     typedef typename E1::matrix_temporary_type result_type;
2896 #endif
2897 };
2898
2899 // (m1 + m2) [i] [j] = m1 [i] [j] + m2 [i] [j]
2900 template<class E1, class E2>
2901 BOOST_UBLAS_INLINE
2902 typename matrix_binary_traits<E1, E2, scalar_plus<typename E1::value_type,
2903 typename E2::value_type> >::result_type
2904 operator + (const matrix_expression<E1> &e1,
2905             const matrix_expression<E2> &e2) {
2906     typedef typename matrix_binary_traits<E1, E2, scalar_plus<typename E1::value_type,
2907         typename E2::value_type> >::expression_type expression_type;
2908     return expression_type (e1 (), e2 ());
2909 }
2910
2911 // (m1 - m2) [i] [j] = m1 [i] [j] - m2 [i] [j]
2912 template<class E1, class E2>
2913 BOOST_UBLAS_INLINE
2914 typename matrix_binary_traits<E1, E2, scalar_minus<typename E1::value_type,
2915 typename E2::value_type> >::result_type
2916 operator - (const matrix_expression<E1> &e1,
2917             const matrix_expression<E2> &e2) {
2918     typedef typename matrix_binary_traits<E1, E2, scalar_minus<typename E1::value_type,
2919         typename E2::value_type> >::expression_type expression_type;
2920     return expression_type (e1 (), e2 ());
2921 }
2922
2923 // (m1 * m2) [i] [j] = m1 [i] [j] * m2 [i] [j]
2924 template<class E1, class E2>
2925 BOOST_UBLAS_INLINE
2926 typename matrix_binary_traits<E1, E2, scalar_multiplies<typename E1::value_type,
2927 typename E2::value_type> >::result_type
2928 element_prod (const matrix_expression<E1> &e1,
2929               const matrix_expression<E2> &e2) {
2930     typedef typename matrix_binary_traits<E1, E2, scalar_multiplies<typename E1::value_type,
2931         typename E2::value_type> >::expression_type expression_type;
2932     return expression_type (e1 (), e2 ());
2933 }
2934
2935 // (m1 / m2) [i] [j] = m1 [i] [j] / m2 [i] [j]
2936 template<class E1, class E2>
2937 BOOST_UBLAS_INLINE
2938 typename matrix_binary_traits<E1, E2, scalar_divides<typename E1::value_type,
2939 typename E2::value_type> >::result_type
2940 element_div (const matrix_expression<E1> &e1,
2941              const matrix_expression<E2> &e2) {
2942     typedef typename matrix_binary_traits<E1, E2, scalar_divides<typename E1::value_type,
2943         typename E2::value_type> >::expression_type expression_type;
2944     return expression_type (e1 (), e2 ());
2945 }
2946
2947 template<class E1, class E2, class F>
2948 class matrix_binary_scalar1:
2949     public matrix_expression<matrix_binary_scalar1<E1, E2, F> > {
2950
2951     typedef E1 expression1_type;
2952     typedef E2 expression2_type;
2953     typedef F functor_type;
2954     typedef const E1& expression1_closure_type;
2955     typedef typename E2::const_closure_type expression2_closure_type;
2956     typedef matrix_binary_scalar1<E1, E2, F> self_type;
2957 public:
2958 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
2959     using matrix_expression<self_type>::operator ();
2960 #endif
2961     typedef typename E2::size_type size_type;
2962     typedef typename E2::difference_type difference_type;
2963     typedef typename F::result_type value_type;
2964     typedef value_type const_reference;
2965     typedef const_reference reference;
2966     typedef const self_type const_closure_type;
2967     typedef const_closure_type closure_type;
2968     typedef typename E2::orientation_category orientation_category;
2969     typedef unknown_storage_tag storage_category;
2970
2971     // Construction and destruction
2972     BOOST_UBLAS_INLINE
2973     matrix_binary_scalar1 (const expression1_type &e1, const expression2_type &e2):
2974       e1_ (e1), e2_ (e2) {}
2975
2976     // Accessors
2977     BOOST_UBLAS_INLINE
2978     size_type size1 () const {
2979         return e2_.size1 ();
2980     }
2981     BOOST_UBLAS_INLINE
2982     size_type size2 () const {
2983         return e2_.size2 ();
2984     }
2985
2986 public:
2987     // Element access
2988     BOOST_UBLAS_INLINE
2989     const_reference operator () (size_type i, size_type j) const {
2990         return functor_type::apply (expression1_type (e1_), e2_ (i, j));
2991     }
2992
2993     // Element access
2994     BOOST_UBLAS_INLINE
2995     const_reference operator () (size_type i) const {
2996         return functor_type::apply (expression1_type (e1_), e2_ (i));
2997     }
2998
2999     // Closure comparison
3000     BOOST_UBLAS_INLINE
3001     bool same_closure (const matrix_binary_scalar1 &mbs1) const {
3002         return &e1_ == &(mbs1.e1_) &&
3003             (*this).e2_.same_closure (mbs1.e2_);
3004     }
3005
3006     // Iterator types
3007 private:
3008     typedef expression1_type const_subiterator1_type;
3009     typedef typename E2::const_iterator1 const_iterator21_type;
3010     typedef typename E2::const_iterator2 const_iterator22_type;
3011     typedef const value_type *const_pointer;
3012
3013 public:
3014 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
3015     typedef indexed_const_iterator1<const_closure_type, typename const_iterator21_type::iterator_category> const_iterator1;
3016     typedef const_iterator1 iterator1;
3017     typedef indexed_const_iterator2<const_closure_type, typename const_iterator22_type::iterator_category> const_iterator2;
3018     typedef const_iterator2 iterator2;
3019 #else
3020     class const_iterator1;
3021     typedef const_iterator1 iterator1;
3022     class const_iterator2;
3023     typedef const_iterator2 iterator2;
3024 #endif
3025     typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
3026     typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
3027
3028     // Element lookup
3029     BOOST_UBLAS_INLINE
3030     const_iterator1 find1 (int rank, size_type i, size_type j) const {
3031         const_iterator21_type it21 (e2_.find1 (rank, i, j));
3032 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
3033         return const_iterator1 (*this, it21.index1 (), it21.index2 ());
3034 #else
3035         return const_iterator1 (*this, const_subiterator1_type (e1_), it21);
3036 #endif
3037     }
3038     BOOST_UBLAS_INLINE
3039     const_iterator2 find2 (int rank, size_type i, size_type j) const {
3040         const_iterator22_type it22 (e2_.find2 (rank, i, j));
3041 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
3042         return const_iterator2 (*this, it22.index1 (), it22.index2 ());
3043 #else
3044         return const_iterator2 (*this, const_subiterator1_type (e1_), it22);
3045 #endif
3046     }
3047
3048     // Iterators enhance the iterators of the referenced expression
3049     // with the binary functor.
3050
3051 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
3052     class const_iterator1:
3053         public container_const_reference<matrix_binary_scalar1>,
3054         public iterator_base_traits<typename E2::const_iterator1::iterator_category>::template
3055         iterator_base<const_iterator1, value_type>::type {
3056     public:
3057         typedef typename E2::const_iterator1::iterator_category iterator_category;
3058         typedef typename matrix_binary_scalar1::difference_type difference_type;
3059         typedef typename matrix_binary_scalar1::value_type value_type;
3060         typedef typename matrix_binary_scalar1::const_reference reference;
3061         typedef typename matrix_binary_scalar1::const_pointer pointer;
3062
3063         typedef const_iterator2 dual_iterator_type;
3064         typedef const_reverse_iterator2 dual_reverse_iterator_type;
3065
3066         // Construction and destruction
3067         BOOST_UBLAS_INLINE
3068         const_iterator1 ():
3069           container_const_reference<self_type> (), it1_ (), it2_ () {}
3070         BOOST_UBLAS_INLINE
3071         const_iterator1 (const self_type &mbs, const const_subiterator1_type &it1, const const_iterator21_type &it2):
3072           container_const_reference<self_type> (mbs), it1_ (it1), it2_ (it2) {}
3073
3074         // Arithmetic
3075         BOOST_UBLAS_INLINE
3076         const_iterator1 &operator ++ () {
3077             ++ it2_;
3078             return *this;
3079         }
3080         BOOST_UBLAS_INLINE
3081         const_iterator1 &operator -- () {
3082             -- it2_ ;
3083             return *this;
3084         }
3085         BOOST_UBLAS_INLINE
3086         const_iterator1 &operator += (difference_type n) {
3087             it2_ += n;
3088             return *this;
3089         }
3090         BOOST_UBLAS_INLINE
3091         const_iterator1 &operator -= (difference_type n) {
3092             it2_ -= n;
3093             return *this;
3094         }
3095         BOOST_UBLAS_INLINE
3096         difference_type operator - (const const_iterator1 &it) const {
3097             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3098             // FIXME we shouldn't compare floats
3099             // BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
3100             return it2_ - it.it2_;
3101         }
3102
3103         // Dereference
3104         BOOST_UBLAS_INLINE
3105         const_reference operator * () const {
3106             return functor_type::apply (it1_, *it2_);
3107         }
3108         BOOST_UBLAS_INLINE
3109         const_reference operator [] (difference_type n) const {
3110             return *(*this + n);
3111         }
3112
3113 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
3114         BOOST_UBLAS_INLINE
3115 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3116         typename self_type::
3117 #endif
3118         const_iterator2 begin () const {
3119             return (*this) ().find2 (1, index1 (), 0);
3120         }
3121         BOOST_UBLAS_INLINE
3122 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3123         typename self_type::
3124 #endif
3125         const_iterator2 cbegin () const {
3126             return begin ();
3127         }
3128         BOOST_UBLAS_INLINE
3129 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3130         typename self_type::
3131 #endif
3132         const_iterator2 end () const {
3133             return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
3134         }
3135         BOOST_UBLAS_INLINE
3136 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3137         typename self_type::
3138 #endif
3139         const_iterator2 cend () const {
3140             return end ();
3141         }
3142         BOOST_UBLAS_INLINE
3143 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3144         typename self_type::
3145 #endif
3146         const_reverse_iterator2 rbegin () const {
3147             return const_reverse_iterator2 (end ());
3148         }
3149         BOOST_UBLAS_INLINE
3150 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3151         typename self_type::
3152 #endif
3153         const_reverse_iterator2 crbegin () const {
3154             return rbegin ();
3155         }
3156         BOOST_UBLAS_INLINE
3157 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3158         typename self_type::
3159 #endif
3160         const_reverse_iterator2 rend () const {
3161             return const_reverse_iterator2 (begin ());
3162         }
3163         BOOST_UBLAS_INLINE
3164 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3165         typename self_type::
3166 #endif
3167         const_reverse_iterator2 crend () const {
3168             return rend ();
3169         }
3170 #endif
3171
3172         // Indices
3173         BOOST_UBLAS_INLINE
3174         size_type index1 () const {
3175             return it2_.index1 ();
3176         }
3177         BOOST_UBLAS_INLINE
3178         size_type index2 () const {
3179             return it2_.index2 ();
3180         }
3181
3182         // Assignment
3183         BOOST_UBLAS_INLINE
3184         const_iterator1 &operator = (const const_iterator1 &it) {
3185             container_const_reference<self_type>::assign (&it ());
3186             it1_ = it.it1_;
3187             it2_ = it.it2_;
3188             return *this;
3189         }
3190
3191         // Comparison
3192         BOOST_UBLAS_INLINE
3193         bool operator == (const const_iterator1 &it) const {
3194             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3195             // FIXME we shouldn't compare floats
3196             // BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
3197             return it2_ == it.it2_;
3198         }
3199         BOOST_UBLAS_INLINE
3200         bool operator < (const const_iterator1 &it) const {
3201             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3202             // FIXME we shouldn't compare floats
3203             // BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
3204             return it2_ < it.it2_;
3205         }
3206
3207     private:
3208         const_subiterator1_type it1_;
3209         const_iterator21_type it2_;
3210     };
3211 #endif
3212
3213     BOOST_UBLAS_INLINE
3214     const_iterator1 begin1 () const {
3215         return find1 (0, 0, 0);
3216     }
3217     BOOST_UBLAS_INLINE
3218     const_iterator1 cbegin1 () const {
3219         return begin1 ();
3220     }
3221     BOOST_UBLAS_INLINE
3222     const_iterator1 end1 () const {
3223         return find1 (0, size1 (), 0);
3224     }
3225     BOOST_UBLAS_INLINE
3226     const_iterator1 cend1 () const {
3227         return end1 ();
3228     }
3229
3230 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
3231     class const_iterator2:
3232         public container_const_reference<matrix_binary_scalar1>,
3233         public iterator_base_traits<typename E2::const_iterator2::iterator_category>::template
3234         iterator_base<const_iterator2, value_type>::type {
3235     public:
3236         typedef typename E2::const_iterator2::iterator_category iterator_category;
3237         typedef typename matrix_binary_scalar1::difference_type difference_type;
3238         typedef typename matrix_binary_scalar1::value_type value_type;
3239         typedef typename matrix_binary_scalar1::const_reference reference;
3240         typedef typename matrix_binary_scalar1::const_pointer pointer;
3241
3242         typedef const_iterator1 dual_iterator_type;
3243         typedef const_reverse_iterator1 dual_reverse_iterator_type;
3244
3245         // Construction and destruction
3246         BOOST_UBLAS_INLINE
3247         const_iterator2 ():
3248           container_const_reference<self_type> (), it1_ (), it2_ () {}
3249         BOOST_UBLAS_INLINE
3250         const_iterator2 (const self_type &mbs, const const_subiterator1_type &it1, const const_iterator22_type &it2):
3251           container_const_reference<self_type> (mbs), it1_ (it1), it2_ (it2) {}
3252
3253         // Arithmetic
3254         BOOST_UBLAS_INLINE
3255         const_iterator2 &operator ++ () {
3256             ++ it2_;
3257             return *this;
3258         }
3259         BOOST_UBLAS_INLINE
3260         const_iterator2 &operator -- () {
3261             -- it2_;
3262             return *this;
3263         }
3264         BOOST_UBLAS_INLINE
3265         const_iterator2 &operator += (difference_type n) {
3266             it2_ += n;
3267             return *this;
3268         }
3269         BOOST_UBLAS_INLINE
3270         const_iterator2 &operator -= (difference_type n) {
3271             it2_ -= n;
3272             return *this;
3273         }
3274         BOOST_UBLAS_INLINE
3275         difference_type operator - (const const_iterator2 &it) const {
3276             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3277             // FIXME we shouldn't compare floats
3278             // BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
3279             return it2_ - it.it2_;
3280         }
3281
3282         // Dereference
3283         BOOST_UBLAS_INLINE
3284         const_reference operator * () const {
3285             return functor_type::apply (it1_, *it2_);
3286         }
3287         BOOST_UBLAS_INLINE
3288         const_reference operator [] (difference_type n) const {
3289             return *(*this + n);
3290         }
3291
3292 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
3293         BOOST_UBLAS_INLINE
3294 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3295         typename self_type::
3296 #endif
3297         const_iterator1 begin () const {
3298             return (*this) ().find1 (1, 0, index2 ());
3299         }
3300         BOOST_UBLAS_INLINE
3301 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3302         typename self_type::
3303 #endif
3304         const_iterator1 cbegin () const {
3305             return begin ();
3306         }
3307         BOOST_UBLAS_INLINE
3308 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3309         typename self_type::
3310 #endif
3311         const_iterator1 end () const {
3312             return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
3313         }
3314         BOOST_UBLAS_INLINE
3315 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3316         typename self_type::
3317 #endif
3318         const_iterator1 cend () const {
3319             return end ();
3320         }
3321         BOOST_UBLAS_INLINE
3322 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3323         typename self_type::
3324 #endif
3325         const_reverse_iterator1 rbegin () const {
3326             return const_reverse_iterator1 (end ());
3327         }
3328         BOOST_UBLAS_INLINE
3329 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3330         typename self_type::
3331 #endif
3332         const_reverse_iterator1 crbegin () const {
3333             return rbegin ();
3334         }
3335         BOOST_UBLAS_INLINE
3336 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3337         typename self_type::
3338 #endif
3339         const_reverse_iterator1 rend () const {
3340             return const_reverse_iterator1 (begin ());
3341         }
3342         BOOST_UBLAS_INLINE
3343 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3344         typename self_type::
3345 #endif
3346         const_reverse_iterator1 crend () const {
3347             return rend ();
3348         }
3349 #endif
3350
3351         // Indices
3352         BOOST_UBLAS_INLINE
3353         size_type index1 () const {
3354             return it2_.index1 ();
3355         }
3356         BOOST_UBLAS_INLINE
3357         size_type index2 () const {
3358             return it2_.index2 ();
3359         }
3360
3361         // Assignment
3362         BOOST_UBLAS_INLINE
3363         const_iterator2 &operator = (const const_iterator2 &it) {
3364             container_const_reference<self_type>::assign (&it ());
3365             it1_ = it.it1_;
3366             it2_ = it.it2_;
3367             return *this;
3368         }
3369
3370         // Comparison
3371         BOOST_UBLAS_INLINE
3372         bool operator == (const const_iterator2 &it) const {
3373             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3374             // FIXME we shouldn't compare floats
3375             // BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
3376             return it2_ == it.it2_;
3377         }
3378         BOOST_UBLAS_INLINE
3379         bool operator < (const const_iterator2 &it) const {
3380             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3381             // FIXME we shouldn't compare floats
3382             // BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
3383             return it2_ < it.it2_;
3384         }
3385
3386     private:
3387         const_subiterator1_type it1_;
3388         const_iterator22_type it2_;
3389     };
3390 #endif
3391
3392     BOOST_UBLAS_INLINE
3393     const_iterator2 begin2 () const {
3394         return find2 (0, 0, 0);
3395     }
3396     BOOST_UBLAS_INLINE
3397     const_iterator2 cbegin2 () const {
3398         return begin2 ();
3399     }
3400     BOOST_UBLAS_INLINE
3401     const_iterator2 end2 () const {
3402         return find2 (0, 0, size2 ());
3403     }
3404     BOOST_UBLAS_INLINE
3405     const_iterator2 cend2 () const {
3406         return end2 ();
3407     }
3408
3409     // Reverse iterators
3410
3411     BOOST_UBLAS_INLINE
3412     const_reverse_iterator1 rbegin1 () const {
3413         return const_reverse_iterator1 (end1 ());
3414     }
3415     BOOST_UBLAS_INLINE
3416     const_reverse_iterator1 crbegin1 () const {
3417         return rbegin1 ();
3418     }
3419     BOOST_UBLAS_INLINE
3420     const_reverse_iterator1 rend1 () const {
3421         return const_reverse_iterator1 (begin1 ());
3422     }
3423     BOOST_UBLAS_INLINE
3424     const_reverse_iterator1 crend1 () const {
3425         return rend1 ();
3426     }
3427
3428     BOOST_UBLAS_INLINE
3429     const_reverse_iterator2 rbegin2 () const {
3430         return const_reverse_iterator2 (end2 ());
3431     }
3432     BOOST_UBLAS_INLINE
3433     const_reverse_iterator2 crbegin2 () const {
3434         return rbegin2 ();
3435     }
3436     BOOST_UBLAS_INLINE
3437     const_reverse_iterator2 rend2 () const {
3438         return const_reverse_iterator2 (begin2 ());
3439     }
3440     BOOST_UBLAS_INLINE
3441     const_reverse_iterator2 crend2 () const {
3442         return rend2 ();
3443     }
3444
3445 private:
3446     expression1_closure_type e1_;
3447     expression2_closure_type e2_;
3448 };
3449
3450 template<class E1, class E2, class F>
3451 struct matrix_binary_scalar1_traits {
3452     typedef matrix_binary_scalar1<E1, E2, F> expression_type;   // allow E1 to be builtin type
3453 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
3454     typedef expression_type result_type;
3455 #else
3456     typedef typename E2::matrix_temporary_type result_type;
3457 #endif
3458 };
3459
3460 // (t * m) [i] [j] = t * m [i] [j]
3461 template<class T1, class E2>
3462 BOOST_UBLAS_INLINE
3463 typename boost::enable_if< is_convertible<T1, typename E2::value_type >,
3464 typename matrix_binary_scalar1_traits<const T1, E2, scalar_multiplies<T1, typename E2::value_type> >::result_type
3465 >::type
3466 operator * (const T1 &e1,
3467             const matrix_expression<E2> &e2) {
3468     typedef typename matrix_binary_scalar1_traits<const T1, E2, scalar_multiplies<T1, typename E2::value_type> >::expression_type expression_type;
3469     return expression_type (e1, e2 ());
3470 }
3471
3472
3473 template<class E1, class E2, class F>
3474 class matrix_binary_scalar2:
3475     public matrix_expression<matrix_binary_scalar2<E1, E2, F> > {
3476
3477     typedef E1 expression1_type;
3478     typedef E2 expression2_type;
3479     typedef F functor_type;
3480 public:
3481     typedef typename E1::const_closure_type expression1_closure_type;
3482     typedef const E2& expression2_closure_type;
3483 private:
3484     typedef matrix_binary_scalar2<E1, E2, F> self_type;
3485 public:
3486 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
3487     using matrix_expression<self_type>::operator ();
3488 #endif
3489     typedef typename E1::size_type size_type;
3490     typedef typename E1::difference_type difference_type;
3491     typedef typename F::result_type value_type;
3492     typedef value_type const_reference;
3493     typedef const_reference reference;
3494
3495     typedef const self_type const_closure_type;
3496     typedef const_closure_type closure_type;
3497     typedef typename E1::orientation_category orientation_category;
3498     typedef unknown_storage_tag storage_category;
3499
3500     // Construction and destruction
3501     BOOST_UBLAS_INLINE
3502     matrix_binary_scalar2 (const expression1_type &e1, const expression2_type &e2):
3503       e1_ (e1), e2_ (e2) {}
3504
3505     // Accessors
3506     BOOST_UBLAS_INLINE
3507     size_type size1 () const {
3508         return e1_.size1 ();
3509     }
3510     BOOST_UBLAS_INLINE
3511     size_type size2 () const {
3512         return e1_.size2 ();
3513     }
3514
3515 public:
3516     // Element access
3517     BOOST_UBLAS_INLINE
3518     const_reference operator () (size_type i, size_type j) const {
3519         return functor_type::apply (e1_ (i, j), expression2_type (e2_));
3520     }
3521     // Element access
3522     BOOST_UBLAS_INLINE
3523     const_reference operator () (size_type i) const {
3524         return functor_type::apply (e1_ (i), expression2_type (e2_));
3525     }
3526
3527     // Closure comparison
3528     BOOST_UBLAS_INLINE
3529     bool same_closure (const matrix_binary_scalar2 &mbs2) const {
3530         return (*this).e1_.same_closure (mbs2.e1_) &&
3531             &e2_ == &(mbs2.e2_);
3532     }
3533
3534     // Iterator types
3535 private:
3536     typedef typename E1::const_iterator1 const_iterator11_type;
3537     typedef typename E1::const_iterator2 const_iterator12_type;
3538     typedef expression2_type const_subiterator2_type;
3539     typedef const value_type *const_pointer;
3540
3541 public:
3542 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
3543     typedef indexed_const_iterator1<const_closure_type, typename const_iterator11_type::iterator_category> const_iterator1;
3544     typedef const_iterator1 iterator1;
3545     typedef indexed_const_iterator2<const_closure_type, typename const_iterator12_type::iterator_category> const_iterator2;
3546     typedef const_iterator2 iterator2;
3547 #else
3548     class const_iterator1;
3549     typedef const_iterator1 iterator1;
3550     class const_iterator2;
3551     typedef const_iterator2 iterator2;
3552 #endif
3553     typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
3554     typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
3555
3556     // Element lookup
3557     BOOST_UBLAS_INLINE
3558     const_iterator1 find1 (int rank, size_type i, size_type j) const {
3559         const_iterator11_type it11 (e1_.find1 (rank, i, j));
3560 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
3561         return const_iterator1 (*this, it11.index1 (), it11.index2 ());
3562 #else
3563         return const_iterator1 (*this, it11, const_subiterator2_type (e2_));
3564 #endif
3565     }
3566     BOOST_UBLAS_INLINE
3567     const_iterator2 find2 (int rank, size_type i, size_type j) const {
3568         const_iterator12_type it12 (e1_.find2 (rank, i, j));
3569 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
3570         return const_iterator2 (*this, it12.index1 (), it12.index2 ());
3571 #else
3572         return const_iterator2 (*this, it12, const_subiterator2_type (e2_));
3573 #endif
3574     }
3575
3576     // Iterators enhance the iterators of the referenced expression
3577     // with the binary functor.
3578
3579 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
3580     class const_iterator1:
3581         public container_const_reference<matrix_binary_scalar2>,
3582         public iterator_base_traits<typename E1::const_iterator1::iterator_category>::template
3583         iterator_base<const_iterator1, value_type>::type {
3584     public:
3585         typedef typename E1::const_iterator1::iterator_category iterator_category;
3586         typedef typename matrix_binary_scalar2::difference_type difference_type;
3587         typedef typename matrix_binary_scalar2::value_type value_type;
3588         typedef typename matrix_binary_scalar2::const_reference reference;
3589         typedef typename matrix_binary_scalar2::const_pointer pointer;
3590
3591         typedef const_iterator2 dual_iterator_type;
3592         typedef const_reverse_iterator2 dual_reverse_iterator_type;
3593
3594         // Construction and destruction
3595         BOOST_UBLAS_INLINE
3596         const_iterator1 ():
3597           container_const_reference<self_type> (), it1_ (), it2_ () {}
3598         BOOST_UBLAS_INLINE
3599         const_iterator1 (const self_type &mbs, const const_iterator11_type &it1, const const_subiterator2_type &it2):
3600           container_const_reference<self_type> (mbs), it1_ (it1), it2_ (it2) {}
3601
3602         // Arithmetic
3603         BOOST_UBLAS_INLINE
3604         const_iterator1 &operator ++ () {
3605             ++ it1_;
3606             return *this;
3607         }
3608         BOOST_UBLAS_INLINE
3609         const_iterator1 &operator -- () {
3610             -- it1_ ;
3611             return *this;
3612         }
3613         BOOST_UBLAS_INLINE
3614         const_iterator1 &operator += (difference_type n) {
3615             it1_ += n;
3616             return *this;
3617         }
3618         BOOST_UBLAS_INLINE
3619         const_iterator1 &operator -= (difference_type n) {
3620             it1_ -= n;
3621             return *this;
3622         }
3623         BOOST_UBLAS_INLINE
3624         difference_type operator - (const const_iterator1 &it) const {
3625             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3626             // FIXME we shouldn't compare floats
3627             // BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
3628             return it1_ - it.it1_;
3629         }
3630
3631         // Dereference
3632         BOOST_UBLAS_INLINE
3633         const_reference operator * () const {
3634             return functor_type::apply (*it1_, it2_);
3635         }
3636         BOOST_UBLAS_INLINE
3637         const_reference operator [] (difference_type n) const {
3638             return *(*this + n);
3639         }
3640
3641 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
3642         BOOST_UBLAS_INLINE
3643 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3644         typename self_type::
3645 #endif
3646         const_iterator2 begin () const {
3647             return (*this) ().find2 (1, index1 (), 0);
3648         }
3649         BOOST_UBLAS_INLINE
3650 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3651         typename self_type::
3652 #endif
3653         const_iterator2 cbegin () const {
3654             return begin ();
3655         }
3656         BOOST_UBLAS_INLINE
3657 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3658         typename self_type::
3659 #endif
3660         const_iterator2 end () const {
3661             return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
3662         }
3663         BOOST_UBLAS_INLINE
3664 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3665         typename self_type::
3666 #endif
3667         const_iterator2 cend () const {
3668             return end ();
3669         }
3670         BOOST_UBLAS_INLINE
3671 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3672         typename self_type::
3673 #endif
3674         const_reverse_iterator2 rbegin () const {
3675             return const_reverse_iterator2 (end ());
3676         }
3677         BOOST_UBLAS_INLINE
3678 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3679         typename self_type::
3680 #endif
3681         const_reverse_iterator2 crbegin () const {
3682             return rbegin ();
3683         }
3684         BOOST_UBLAS_INLINE
3685 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3686         typename self_type::
3687 #endif
3688         const_reverse_iterator2 rend () const {
3689             return const_reverse_iterator2 (begin ());
3690         }
3691         BOOST_UBLAS_INLINE
3692 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3693         typename self_type::
3694 #endif
3695         const_reverse_iterator2 crend () const {
3696             return rend ();
3697         }
3698 #endif
3699
3700         // Indices
3701         BOOST_UBLAS_INLINE
3702         size_type index1 () const {
3703             return it1_.index1 ();
3704         }
3705         BOOST_UBLAS_INLINE
3706         size_type index2 () const {
3707             return it1_.index2 ();
3708         }
3709
3710         // Assignment
3711         BOOST_UBLAS_INLINE
3712         const_iterator1 &operator = (const const_iterator1 &it) {
3713             container_const_reference<self_type>::assign (&it ());
3714             it1_ = it.it1_;
3715             it2_ = it.it2_;
3716             return *this;
3717         }
3718
3719         // Comparison
3720         BOOST_UBLAS_INLINE
3721         bool operator == (const const_iterator1 &it) const {
3722             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3723             // FIXME we shouldn't compare floats
3724             // BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
3725             return it1_ == it.it1_;
3726         }
3727         BOOST_UBLAS_INLINE
3728         bool operator < (const const_iterator1 &it) const {
3729             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3730             // FIXME we shouldn't compare floats
3731             // BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
3732             return it1_ < it.it1_;
3733         }
3734
3735     private:
3736         const_iterator11_type it1_;
3737         const_subiterator2_type it2_;
3738     };
3739 #endif
3740
3741     BOOST_UBLAS_INLINE
3742     const_iterator1 begin1 () const {
3743         return find1 (0, 0, 0);
3744     }
3745     BOOST_UBLAS_INLINE
3746     const_iterator1 cbegin1 () const {
3747         return begin1 ();
3748     }
3749     BOOST_UBLAS_INLINE
3750     const_iterator1 end1 () const {
3751         return find1 (0, size1 (), 0);
3752     }
3753     BOOST_UBLAS_INLINE
3754     const_iterator1 cend1 () const {
3755         return end1 ();
3756     }
3757
3758 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
3759     class const_iterator2:
3760         public container_const_reference<matrix_binary_scalar2>,
3761         public iterator_base_traits<typename E1::const_iterator2::iterator_category>::template
3762         iterator_base<const_iterator2, value_type>::type {
3763     public:
3764         typedef typename E1::const_iterator2::iterator_category iterator_category;
3765         typedef typename matrix_binary_scalar2::difference_type difference_type;
3766         typedef typename matrix_binary_scalar2::value_type value_type;
3767         typedef typename matrix_binary_scalar2::const_reference reference;
3768         typedef typename matrix_binary_scalar2::const_pointer pointer;
3769
3770         typedef const_iterator1 dual_iterator_type;
3771         typedef const_reverse_iterator1 dual_reverse_iterator_type;
3772
3773         // Construction and destruction
3774         BOOST_UBLAS_INLINE
3775         const_iterator2 ():
3776           container_const_reference<self_type> (), it1_ (), it2_ () {}
3777         BOOST_UBLAS_INLINE
3778         const_iterator2 (const self_type &mbs, const const_iterator12_type &it1, const const_subiterator2_type &it2):
3779           container_const_reference<self_type> (mbs), it1_ (it1), it2_ (it2) {}
3780
3781         // Arithmetic
3782         BOOST_UBLAS_INLINE
3783         const_iterator2 &operator ++ () {
3784             ++ it1_;
3785             return *this;
3786         }
3787         BOOST_UBLAS_INLINE
3788         const_iterator2 &operator -- () {
3789             -- it1_;
3790             return *this;
3791         }
3792         BOOST_UBLAS_INLINE
3793         const_iterator2 &operator += (difference_type n) {
3794             it1_ += n;
3795             return *this;
3796         }
3797         BOOST_UBLAS_INLINE
3798         const_iterator2 &operator -= (difference_type n) {
3799             it1_ -= n;
3800             return *this;
3801         }
3802         BOOST_UBLAS_INLINE
3803         difference_type operator - (const const_iterator2 &it) const {
3804             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3805             // FIXME we shouldn't compare floats
3806             // BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
3807             return it1_ - it.it1_;
3808         }
3809
3810         // Dereference
3811         BOOST_UBLAS_INLINE
3812         const_reference operator * () const {
3813             return functor_type::apply (*it1_, it2_);
3814         }
3815         BOOST_UBLAS_INLINE
3816         const_reference operator [] (difference_type n) const {
3817             return *(*this + n);
3818         }
3819
3820 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
3821         BOOST_UBLAS_INLINE
3822 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3823         typename self_type::
3824 #endif
3825         const_iterator1 begin () const {
3826             return (*this) ().find1 (1, 0, index2 ());
3827         }
3828         BOOST_UBLAS_INLINE
3829 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3830         typename self_type::
3831 #endif
3832         const_iterator1 cbegin () const {
3833             return begin ();
3834         }
3835         BOOST_UBLAS_INLINE
3836 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3837         typename self_type::
3838 #endif
3839         const_iterator1 end () const {
3840             return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
3841         }
3842         BOOST_UBLAS_INLINE
3843 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3844         typename self_type::
3845 #endif
3846         const_iterator1 cend () const {
3847             return end ();
3848         }
3849         BOOST_UBLAS_INLINE
3850 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3851         typename self_type::
3852 #endif
3853         const_reverse_iterator1 rbegin () const {
3854             return const_reverse_iterator1 (end ());
3855         }
3856         BOOST_UBLAS_INLINE
3857 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3858         typename self_type::
3859 #endif
3860         const_reverse_iterator1 crbegin () const {
3861             return rbegin ();
3862         }
3863         BOOST_UBLAS_INLINE
3864 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3865         typename self_type::
3866 #endif
3867         const_reverse_iterator1 rend () const {
3868             return const_reverse_iterator1 (begin ());
3869         }
3870         BOOST_UBLAS_INLINE
3871 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3872         typename self_type::
3873 #endif
3874         const_reverse_iterator1 crend () const {
3875             return rend ();
3876         }
3877 #endif
3878
3879         // Indices
3880         BOOST_UBLAS_INLINE
3881         size_type index1 () const {
3882             return it1_.index1 ();
3883         }
3884         BOOST_UBLAS_INLINE
3885         size_type index2 () const {
3886             return it1_.index2 ();
3887         }
3888
3889         // Assignment
3890         BOOST_UBLAS_INLINE
3891         const_iterator2 &operator = (const const_iterator2 &it) {
3892             container_const_reference<self_type>::assign (&it ());
3893             it1_ = it.it1_;
3894             it2_ = it.it2_;
3895             return *this;
3896         }
3897
3898         // Comparison
3899         BOOST_UBLAS_INLINE
3900         bool operator == (const const_iterator2 &it) const {
3901             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3902             // FIXME we shouldn't compare floats
3903             // BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
3904             return it1_ == it.it1_;
3905         }
3906         BOOST_UBLAS_INLINE
3907         bool operator < (const const_iterator2 &it) const {
3908             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3909             // FIXME we shouldn't compare floats
3910             // BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
3911             return it1_ < it.it1_;
3912         }
3913
3914     private:
3915         const_iterator12_type it1_;
3916         const_subiterator2_type it2_;
3917     };
3918 #endif
3919
3920     BOOST_UBLAS_INLINE
3921     const_iterator2 begin2 () const {
3922         return find2 (0, 0, 0);
3923     }
3924     BOOST_UBLAS_INLINE
3925     const_iterator2 cbegin2 () const {
3926         return begin2 ();
3927     }
3928     BOOST_UBLAS_INLINE
3929     const_iterator2 end2 () const {
3930         return find2 (0, 0, size2 ());
3931     }
3932     BOOST_UBLAS_INLINE
3933     const_iterator2 cend2 () const {
3934         return end2 ();
3935     }
3936
3937     // Reverse iterators
3938
3939     BOOST_UBLAS_INLINE
3940     const_reverse_iterator1 rbegin1 () const {
3941         return const_reverse_iterator1 (end1 ());
3942     }
3943     BOOST_UBLAS_INLINE
3944     const_reverse_iterator1 crbegin1 () const {
3945         return rbegin1 ();
3946     }
3947     BOOST_UBLAS_INLINE
3948     const_reverse_iterator1 rend1 () const {
3949         return const_reverse_iterator1 (begin1 ());
3950     }
3951     BOOST_UBLAS_INLINE
3952     const_reverse_iterator1 crend1 () const {
3953         return rend1 ();
3954     }
3955
3956     BOOST_UBLAS_INLINE
3957     const_reverse_iterator2 rbegin2 () const {
3958         return const_reverse_iterator2 (end2 ());
3959     }
3960     BOOST_UBLAS_INLINE
3961     const_reverse_iterator2 crbegin2 () const {
3962         return rbegin2 ();
3963     }
3964     BOOST_UBLAS_INLINE
3965     const_reverse_iterator2 rend2 () const {
3966         return const_reverse_iterator2 (begin2 ());
3967     }
3968     BOOST_UBLAS_INLINE
3969     const_reverse_iterator2 crend2 () const {
3970         return rend2 ();
3971     }
3972
3973 private:
3974     expression1_closure_type e1_;
3975     expression2_closure_type e2_;
3976 };
3977
3978 template<class E1, class E2, class F>
3979 struct matrix_binary_scalar2_traits {
3980     typedef matrix_binary_scalar2<E1, E2, F> expression_type;   // allow E2 to be builtin type
3981 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
3982     typedef expression_type result_type;
3983 #else
3984     typedef typename E1::matrix_temporary_type result_type;
3985 #endif
3986 };
3987
3988 // (m * t) [i] [j] = m [i] [j] * t
3989 template<class E1, class T2>
3990 BOOST_UBLAS_INLINE
3991 typename boost::enable_if< is_convertible<T2, typename E1::value_type>,
3992 typename matrix_binary_scalar2_traits<E1, const T2, scalar_multiplies<typename E1::value_type, T2> >::result_type
3993 >::type
3994 operator * (const matrix_expression<E1> &e1,
3995             const T2 &e2) {
3996     typedef typename matrix_binary_scalar2_traits<E1, const T2, scalar_multiplies<typename E1::value_type, T2> >::expression_type expression_type;
3997     return expression_type (e1 (), e2);
3998 }
3999
4000 // (m / t) [i] [j] = m [i] [j] / t
4001 template<class E1, class T2>
4002 BOOST_UBLAS_INLINE
4003 typename boost::enable_if< is_convertible<T2, typename E1::value_type>,
4004 typename matrix_binary_scalar2_traits<E1, const T2, scalar_divides<typename E1::value_type, T2> >::result_type
4005 >::type
4006 operator / (const matrix_expression<E1> &e1,
4007             const T2 &e2) {
4008     typedef typename matrix_binary_scalar2_traits<E1, const T2, scalar_divides<typename E1::value_type, T2> >::expression_type expression_type;
4009     return expression_type (e1 (), e2);
4010 }
4011
4012
4013 template<class E1, class E2, class F>
4014 class matrix_vector_binary1:
4015     public vector_expression<matrix_vector_binary1<E1, E2, F> > {
4016
4017 public:
4018     typedef E1 expression1_type;
4019     typedef E2 expression2_type;
4020 private:
4021     typedef F functor_type;
4022 public:
4023     typedef typename E1::const_closure_type expression1_closure_type;
4024     typedef typename E2::const_closure_type expression2_closure_type;
4025 private:
4026     typedef matrix_vector_binary1<E1, E2, F> self_type;
4027 public:
4028 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
4029     using vector_expression<self_type>::operator ();
4030 #endif
4031     static const unsigned complexity = 1;
4032     typedef typename promote_traits<typename E1::size_type, typename E2::size_type>::promote_type size_type;
4033     typedef typename promote_traits<typename E1::difference_type, typename E2::difference_type>::promote_type difference_type;
4034     typedef typename F::result_type value_type;
4035     typedef value_type const_reference;
4036     typedef const_reference reference;
4037     typedef const self_type const_closure_type;
4038     typedef const_closure_type closure_type;
4039     typedef unknown_storage_tag storage_category;
4040
4041     // Construction and destruction
4042     BOOST_UBLAS_INLINE
4043     matrix_vector_binary1 (const expression1_type &e1, const expression2_type &e2):
4044       e1_ (e1), e2_ (e2) {}
4045
4046     // Accessors
4047     BOOST_UBLAS_INLINE
4048     size_type size () const {
4049         return e1_.size1 ();
4050     }
4051
4052 public:
4053     // Expression accessors
4054     BOOST_UBLAS_INLINE
4055     const expression1_closure_type &expression1 () const {
4056         return e1_;
4057     }
4058     BOOST_UBLAS_INLINE
4059     const expression2_closure_type &expression2 () const {
4060         return e2_;
4061     }
4062
4063 public:
4064     // Element access
4065     BOOST_UBLAS_INLINE
4066     const_reference operator () (size_type i) const {
4067         return functor_type::apply (e1_, e2_, i);
4068     }
4069
4070     // Closure comparison
4071     BOOST_UBLAS_INLINE
4072     bool same_closure (const matrix_vector_binary1 &mvb1) const {
4073         return (*this).expression1 ().same_closure (mvb1.expression1 ()) &&
4074             (*this).expression2 ().same_closure (mvb1.expression2 ());
4075     }
4076
4077     // Iterator types
4078 private:
4079     typedef typename E1::const_iterator1 const_subiterator1_type;
4080     typedef typename E2::const_iterator const_subiterator2_type;
4081     typedef const value_type *const_pointer;
4082
4083 public:
4084 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
4085     typedef indexed_const_iterator<const_closure_type, typename const_subiterator1_type::iterator_category> const_iterator;
4086     typedef const_iterator iterator;
4087 #else
4088     class const_iterator;
4089     typedef const_iterator iterator;
4090 #endif
4091
4092     // Element lookup
4093     BOOST_UBLAS_INLINE
4094     const_iterator find (size_type i) const {
4095 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
4096         const_subiterator1_type it1 (e1_.find1 (0, i, 0));
4097         return const_iterator (*this, it1.index1 ());
4098 #else
4099         return const_iterator (*this, e1_.find1 (0, i, 0));
4100 #endif
4101     }
4102
4103
4104 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
4105     class const_iterator:
4106         public container_const_reference<matrix_vector_binary1>,
4107         public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
4108         typename E2::const_iterator::iterator_category>::iterator_category>::template
4109         iterator_base<const_iterator, value_type>::type {
4110     public:
4111         typedef typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
4112         typename E2::const_iterator::iterator_category>::iterator_category iterator_category;
4113         typedef typename matrix_vector_binary1::difference_type difference_type;
4114         typedef typename matrix_vector_binary1::value_type value_type;
4115         typedef typename matrix_vector_binary1::const_reference reference;
4116         typedef typename matrix_vector_binary1::const_pointer pointer;
4117
4118         // Construction and destruction
4119 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4120         BOOST_UBLAS_INLINE
4121         const_iterator ():
4122           container_const_reference<self_type> (), it1_ (), e2_begin_ (), e2_end_ () {}
4123         BOOST_UBLAS_INLINE
4124         const_iterator (const self_type &mvb, const const_subiterator1_type &it1):
4125           container_const_reference<self_type> (mvb), it1_ (it1), e2_begin_ (mvb.expression2 ().begin ()), e2_end_ (mvb.expression2 ().end ()) {}
4126 #else
4127         BOOST_UBLAS_INLINE
4128         const_iterator ():
4129           container_const_reference<self_type> (), it1_ () {}
4130         BOOST_UBLAS_INLINE
4131         const_iterator (const self_type &mvb, const const_subiterator1_type &it1):
4132           container_const_reference<self_type> (mvb), it1_ (it1) {}
4133 #endif
4134
4135     private:
4136         // Dense random access specialization
4137         BOOST_UBLAS_INLINE
4138         value_type dereference (dense_random_access_iterator_tag) const {
4139             const self_type &mvb = (*this) ();
4140 #ifdef BOOST_UBLAS_USE_INDEXING
4141             return mvb (index ());
4142 #elif BOOST_UBLAS_USE_ITERATING
4143             difference_type size = BOOST_UBLAS_SAME (mvb.expression1 ().size2 (), mvb.expression2 ().size ());
4144 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4145             return functor_type::apply (size, it1_.begin (), e2_begin_);
4146 #else
4147             return functor_type::apply (size, it1_.begin (), mvb.expression2 ().begin ());
4148 #endif
4149 #else
4150             difference_type size = BOOST_UBLAS_SAME (mvb.expression1 ().size2 (), mvb.expression2 ().size ());
4151             if (size >= BOOST_UBLAS_ITERATOR_THRESHOLD)
4152 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4153                 return functor_type::apply (size, it1_.begin (), e2_begin_);
4154 #else
4155                 return functor_type::apply (size, it1_.begin (), mvb.expression2 ().begin ());
4156 #endif
4157             else
4158                 return mvb (index ());
4159 #endif
4160         }
4161
4162         // Packed bidirectional specialization
4163         BOOST_UBLAS_INLINE
4164         value_type dereference (packed_random_access_iterator_tag) const {
4165 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4166             return functor_type::apply (it1_.begin (), it1_.end (), e2_begin_, e2_end_);
4167 #else
4168             const self_type &mvb = (*this) ();
4169 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
4170             return functor_type::apply (it1_.begin (), it1_.end (),
4171                                         mvb.expression2 ().begin (), mvb.expression2 ().end ());
4172 #else
4173             return functor_type::apply (boost::numeric::ublas::begin (it1_, iterator1_tag ()),
4174                                         boost::numeric::ublas::end (it1_, iterator1_tag ()),
4175                                         mvb.expression2 ().begin (), mvb.expression2 ().end ());
4176 #endif
4177 #endif
4178         }
4179
4180         // Sparse bidirectional specialization
4181         BOOST_UBLAS_INLINE
4182         value_type dereference (sparse_bidirectional_iterator_tag) const {
4183 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4184             return functor_type::apply (it1_.begin (), it1_.end (), e2_begin_, e2_end_, sparse_bidirectional_iterator_tag ());
4185 #else
4186             const self_type &mvb = (*this) ();
4187 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
4188             return functor_type::apply (it1_.begin (), it1_.end (),
4189                                         mvb.expression2 ().begin (), mvb.expression2 ().end (), sparse_bidirectional_iterator_tag ());
4190 #else
4191             return functor_type::apply (boost::numeric::ublas::begin (it1_, iterator1_tag ()),
4192                                         boost::numeric::ublas::end (it1_, iterator1_tag ()),
4193                                         mvb.expression2 ().begin (), mvb.expression2 ().end (), sparse_bidirectional_iterator_tag ());
4194 #endif
4195 #endif
4196         }
4197
4198     public:
4199         // Arithmetic
4200         BOOST_UBLAS_INLINE
4201         const_iterator &operator ++ () {
4202             ++ it1_;
4203             return *this;
4204         }
4205         BOOST_UBLAS_INLINE
4206         const_iterator &operator -- () {
4207             -- it1_;
4208             return *this;
4209         }
4210         BOOST_UBLAS_INLINE
4211         const_iterator &operator += (difference_type n) {
4212             it1_ += n;
4213             return *this;
4214         }
4215         BOOST_UBLAS_INLINE
4216         const_iterator &operator -= (difference_type n) {
4217             it1_ -= n;
4218             return *this;
4219         }
4220         BOOST_UBLAS_INLINE
4221         difference_type operator - (const const_iterator &it) const {
4222             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
4223             return it1_ - it.it1_;
4224         }
4225
4226         // Dereference
4227         BOOST_UBLAS_INLINE
4228         const_reference operator * () const {
4229             return dereference (iterator_category ());
4230         }
4231         BOOST_UBLAS_INLINE
4232         const_reference operator [] (difference_type n) const {
4233             return *(*this + n);
4234         }
4235
4236         // Index
4237         BOOST_UBLAS_INLINE
4238         size_type index () const {
4239             return it1_.index1 ();
4240         }
4241
4242         // Assignment
4243         BOOST_UBLAS_INLINE
4244         const_iterator &operator = (const const_iterator &it) {
4245             container_const_reference<self_type>::assign (&it ());
4246             it1_ = it.it1_;
4247 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4248             e2_begin_ = it.e2_begin_;
4249             e2_end_ = it.e2_end_;
4250 #endif
4251             return *this;
4252         }
4253
4254         // Comparison
4255         BOOST_UBLAS_INLINE
4256         bool operator == (const const_iterator &it) const {
4257             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
4258             return it1_ == it.it1_;
4259         }
4260         BOOST_UBLAS_INLINE
4261         bool operator < (const const_iterator &it) const {
4262             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
4263             return it1_ < it.it1_;
4264         }
4265
4266     private:
4267         const_subiterator1_type it1_;
4268 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4269         // Mutable due to assignment
4270         /* const */ const_subiterator2_type e2_begin_;
4271         /* const */ const_subiterator2_type e2_end_;
4272 #endif
4273     };
4274 #endif
4275
4276     BOOST_UBLAS_INLINE
4277     const_iterator begin () const {
4278         return find (0);
4279     }
4280     BOOST_UBLAS_INLINE
4281     const_iterator cbegin () const {
4282         return begin ();
4283     }
4284     BOOST_UBLAS_INLINE
4285     const_iterator end () const {
4286         return find (size ());
4287     }
4288     BOOST_UBLAS_INLINE
4289     const_iterator cend () const {
4290         return end ();
4291     }
4292
4293     // Reverse iterator
4294     typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
4295
4296     BOOST_UBLAS_INLINE
4297     const_reverse_iterator rbegin () const {
4298         return const_reverse_iterator (end ());
4299     }
4300     BOOST_UBLAS_INLINE
4301     const_reverse_iterator crbegin () const {
4302         return rbegin ();
4303     }
4304     BOOST_UBLAS_INLINE
4305     const_reverse_iterator rend () const {
4306         return const_reverse_iterator (begin ());
4307     }
4308     BOOST_UBLAS_INLINE
4309     const_reverse_iterator crend () const {
4310         return rend ();
4311     }
4312
4313 private:
4314     expression1_closure_type e1_;
4315     expression2_closure_type e2_;
4316 };
4317
4318 template<class T1, class E1, class T2, class E2>
4319 struct matrix_vector_binary1_traits {
4320     typedef unknown_storage_tag storage_category;
4321     typedef row_major_tag orientation_category;
4322     typedef typename promote_traits<T1, T2>::promote_type promote_type;
4323     typedef matrix_vector_binary1<E1, E2, matrix_vector_prod1<E1, E2, promote_type> > expression_type;
4324 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
4325     typedef expression_type result_type;
4326 #else
4327     typedef typename E1::vector_temporary_type result_type;
4328 #endif
4329 };
4330
4331 template<class E1, class E2>
4332 BOOST_UBLAS_INLINE
4333 typename matrix_vector_binary1_traits<typename E1::value_type, E1,
4334 typename E2::value_type, E2>::result_type
4335 prod (const matrix_expression<E1> &e1,
4336       const vector_expression<E2> &e2,
4337       unknown_storage_tag,
4338       row_major_tag) {
4339     typedef typename matrix_vector_binary1_traits<typename E1::value_type, E1,
4340         typename E2::value_type, E2>::expression_type expression_type;
4341     return expression_type (e1 (), e2 ());
4342 }
4343
4344 // Dispatcher
4345 template<class E1, class E2>
4346 BOOST_UBLAS_INLINE
4347 typename matrix_vector_binary1_traits<typename E1::value_type, E1,
4348 typename E2::value_type, E2>::result_type
4349 prod (const matrix_expression<E1> &e1,
4350       const vector_expression<E2> &e2) {
4351     BOOST_STATIC_ASSERT (E2::complexity == 0);
4352     typedef typename matrix_vector_binary1_traits<typename E1::value_type, E1,
4353         typename E2::value_type, E2>::storage_category storage_category;
4354     typedef typename matrix_vector_binary1_traits<typename E1::value_type, E1,
4355         typename E2::value_type, E2>::orientation_category orientation_category;
4356     return prod (e1, e2, storage_category (), orientation_category ());
4357 }
4358
4359 template<class E1, class E2>
4360 BOOST_UBLAS_INLINE
4361 typename matrix_vector_binary1_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4362 typename type_traits<typename E2::value_type>::precision_type, E2>::result_type
4363 prec_prod (const matrix_expression<E1> &e1,
4364            const vector_expression<E2> &e2,
4365            unknown_storage_tag,
4366            row_major_tag) {
4367     typedef typename matrix_vector_binary1_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4368         typename type_traits<typename E2::value_type>::precision_type, E2>::expression_type expression_type;
4369     return expression_type (e1 (), e2 ());
4370 }
4371
4372 // Dispatcher
4373 template<class E1, class E2>
4374 BOOST_UBLAS_INLINE
4375 typename matrix_vector_binary1_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4376 typename type_traits<typename E2::value_type>::precision_type, E2>::result_type
4377 prec_prod (const matrix_expression<E1> &e1,
4378            const vector_expression<E2> &e2) {
4379     BOOST_STATIC_ASSERT (E2::complexity == 0);
4380     typedef typename matrix_vector_binary1_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4381         typename type_traits<typename E2::value_type>::precision_type, E2>::storage_category storage_category;
4382     typedef typename matrix_vector_binary1_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4383         typename type_traits<typename E2::value_type>::precision_type, E2>::orientation_category orientation_category;
4384     return prec_prod (e1, e2, storage_category (), orientation_category ());
4385 }
4386
4387 template<class V, class E1, class E2>
4388 BOOST_UBLAS_INLINE
4389 V &
4390 prod (const matrix_expression<E1> &e1,
4391       const vector_expression<E2> &e2,
4392       V &v) {
4393     return v.assign (prod (e1, e2));
4394 }
4395
4396 template<class V, class E1, class E2>
4397 BOOST_UBLAS_INLINE
4398 V &
4399 prec_prod (const matrix_expression<E1> &e1,
4400            const vector_expression<E2> &e2,
4401            V &v) {
4402     return v.assign (prec_prod (e1, e2));
4403 }
4404
4405 template<class V, class E1, class E2>
4406 BOOST_UBLAS_INLINE
4407 V
4408 prod (const matrix_expression<E1> &e1,
4409       const vector_expression<E2> &e2) {
4410     return V (prod (e1, e2));
4411 }
4412
4413 template<class V, class E1, class E2>
4414 BOOST_UBLAS_INLINE
4415 V
4416 prec_prod (const matrix_expression<E1> &e1,
4417            const vector_expression<E2> &e2) {
4418     return V (prec_prod (e1, e2));
4419 }
4420
4421 template<class E1, class E2, class F>
4422 class matrix_vector_binary2:
4423     public vector_expression<matrix_vector_binary2<E1, E2, F> > {
4424
4425     typedef E1 expression1_type;
4426     typedef E2 expression2_type;
4427     typedef F functor_type;
4428 public:
4429     typedef typename E1::const_closure_type expression1_closure_type;
4430     typedef typename E2::const_closure_type expression2_closure_type;
4431 private:
4432     typedef matrix_vector_binary2<E1, E2, F> self_type;
4433 public:
4434 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
4435     using vector_expression<self_type>::operator ();
4436 #endif
4437     static const unsigned complexity = 1;
4438     typedef typename promote_traits<typename E1::size_type, typename E2::size_type>::promote_type size_type;
4439     typedef typename promote_traits<typename E1::difference_type, typename E2::difference_type>::promote_type difference_type;
4440     typedef typename F::result_type value_type;
4441     typedef value_type const_reference;
4442     typedef const_reference reference;
4443     typedef const self_type const_closure_type;
4444     typedef const_closure_type closure_type;
4445     typedef unknown_storage_tag storage_category;
4446
4447     // Construction and destruction
4448     BOOST_UBLAS_INLINE
4449     matrix_vector_binary2 (const expression1_type &e1, const expression2_type &e2):
4450       e1_ (e1), e2_ (e2) {}
4451
4452     // Accessors
4453     BOOST_UBLAS_INLINE
4454     size_type size () const {
4455         return e2_.size2 ();
4456     }
4457
4458 public:
4459     // Expression accessors
4460     BOOST_UBLAS_INLINE
4461     const expression1_closure_type &expression1 () const {
4462         return e1_;
4463     }
4464     BOOST_UBLAS_INLINE
4465     const expression2_closure_type &expression2 () const {
4466         return e2_;
4467     }
4468 public:
4469
4470     // Element access
4471     BOOST_UBLAS_INLINE
4472     const_reference operator () (size_type j) const {
4473         return functor_type::apply (e1_, e2_, j);
4474     }
4475
4476     // Closure comparison
4477     BOOST_UBLAS_INLINE
4478     bool same_closure (const matrix_vector_binary2 &mvb2) const {
4479         return (*this).expression1 ().same_closure (mvb2.expression1 ()) &&
4480             (*this).expression2 ().same_closure (mvb2.expression2 ());
4481     }
4482
4483     // Iterator types
4484 private:
4485     typedef typename E1::const_iterator const_subiterator1_type;
4486     typedef typename E2::const_iterator2 const_subiterator2_type;
4487     typedef const value_type *const_pointer;
4488
4489 public:
4490 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
4491     typedef indexed_const_iterator<const_closure_type, typename const_subiterator2_type::iterator_category> const_iterator;
4492     typedef const_iterator iterator;
4493 #else
4494     class const_iterator;
4495     typedef const_iterator iterator;
4496 #endif
4497
4498     // Element lookup
4499     BOOST_UBLAS_INLINE
4500     const_iterator find (size_type j) const {
4501 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
4502         const_subiterator2_type it2 (e2_.find2 (0, 0, j));
4503         return const_iterator (*this, it2.index2 ());
4504 #else
4505         return const_iterator (*this, e2_.find2 (0, 0, j));
4506 #endif
4507     }
4508
4509
4510 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
4511     class const_iterator:
4512         public container_const_reference<matrix_vector_binary2>,
4513         public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
4514         typename E2::const_iterator2::iterator_category>::iterator_category>::template
4515         iterator_base<const_iterator, value_type>::type {
4516     public:
4517         typedef typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
4518         typename E2::const_iterator2::iterator_category>::iterator_category iterator_category;
4519         typedef typename matrix_vector_binary2::difference_type difference_type;
4520         typedef typename matrix_vector_binary2::value_type value_type;
4521         typedef typename matrix_vector_binary2::const_reference reference;
4522         typedef typename matrix_vector_binary2::const_pointer pointer;
4523
4524         // Construction and destruction
4525 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4526         BOOST_UBLAS_INLINE
4527         const_iterator ():
4528           container_const_reference<self_type> (), it2_ (), e1_begin_ (), e1_end_ () {}
4529         BOOST_UBLAS_INLINE
4530         const_iterator (const self_type &mvb, const const_subiterator2_type &it2):
4531           container_const_reference<self_type> (mvb), it2_ (it2), e1_begin_ (mvb.expression1 ().begin ()), e1_end_ (mvb.expression1 ().end ()) {}
4532 #else
4533         BOOST_UBLAS_INLINE
4534         const_iterator ():
4535           container_const_reference<self_type> (), it2_ () {}
4536         BOOST_UBLAS_INLINE
4537         const_iterator (const self_type &mvb, const const_subiterator2_type &it2):
4538           container_const_reference<self_type> (mvb), it2_ (it2) {}
4539 #endif
4540
4541     private:
4542         // Dense random access specialization
4543         BOOST_UBLAS_INLINE
4544         value_type dereference (dense_random_access_iterator_tag) const {
4545             const self_type &mvb = (*this) ();
4546 #ifdef BOOST_UBLAS_USE_INDEXING
4547             return mvb (index ());
4548 #elif BOOST_UBLAS_USE_ITERATING
4549             difference_type size = BOOST_UBLAS_SAME (mvb.expression2 ().size1 (), mvb.expression1 ().size ());
4550 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4551             return functor_type::apply (size, e1_begin_, it2_.begin ());
4552 #else
4553             return functor_type::apply (size, mvb.expression1 ().begin (), it2_.begin ());
4554 #endif
4555 #else
4556             difference_type size = BOOST_UBLAS_SAME (mvb.expression2 ().size1 (), mvb.expression1 ().size ());
4557             if (size >= BOOST_UBLAS_ITERATOR_THRESHOLD)
4558 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4559                 return functor_type::apply (size, e1_begin_, it2_.begin ());
4560 #else
4561                 return functor_type::apply (size, mvb.expression1 ().begin (), it2_.begin ());
4562 #endif
4563             else
4564                 return mvb (index ());
4565 #endif
4566         }
4567
4568         // Packed bidirectional specialization
4569         BOOST_UBLAS_INLINE
4570         value_type dereference (packed_random_access_iterator_tag) const {
4571 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4572             return functor_type::apply (e1_begin_, e1_end_, it2_.begin (), it2_.end ());
4573 #else
4574             const self_type &mvb = (*this) ();
4575 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
4576             return functor_type::apply (mvb.expression1 ().begin (), mvb.expression1 ().end (),
4577                                         it2_.begin (), it2_.end ());
4578 #else
4579             return functor_type::apply (mvb.expression1 ().begin (), mvb.expression1 ().end (),
4580                                         boost::numeric::ublas::begin (it2_, iterator2_tag ()),
4581                                         boost::numeric::ublas::end (it2_, iterator2_tag ()));
4582 #endif
4583 #endif
4584         }
4585
4586         // Sparse bidirectional specialization
4587         BOOST_UBLAS_INLINE
4588         value_type dereference (sparse_bidirectional_iterator_tag) const {
4589 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4590             return functor_type::apply (e1_begin_, e1_end_, it2_.begin (), it2_.end (), sparse_bidirectional_iterator_tag ());
4591 #else
4592             const self_type &mvb = (*this) ();
4593 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
4594             return functor_type::apply (mvb.expression1 ().begin (), mvb.expression1 ().end (),
4595                                         it2_.begin (), it2_.end (), sparse_bidirectional_iterator_tag ());
4596 #else
4597             return functor_type::apply (mvb.expression1 ().begin (), mvb.expression1 ().end (),
4598                                         boost::numeric::ublas::begin (it2_, iterator2_tag ()),
4599                                         boost::numeric::ublas::end (it2_, iterator2_tag ()), sparse_bidirectional_iterator_tag ());
4600 #endif
4601 #endif
4602         }
4603
4604     public:
4605         // Arithmetic
4606         BOOST_UBLAS_INLINE
4607         const_iterator &operator ++ () {
4608             ++ it2_;
4609             return *this;
4610         }
4611         BOOST_UBLAS_INLINE
4612         const_iterator &operator -- () {
4613             -- it2_;
4614             return *this;
4615         }
4616         BOOST_UBLAS_INLINE
4617         const_iterator &operator += (difference_type n) {
4618             it2_ += n;
4619             return *this;
4620         }
4621         BOOST_UBLAS_INLINE
4622         const_iterator &operator -= (difference_type n) {
4623             it2_ -= n;
4624             return *this;
4625         }
4626         BOOST_UBLAS_INLINE
4627         difference_type operator - (const const_iterator &it) const {
4628             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
4629             return it2_ - it.it2_;
4630         }
4631
4632         // Dereference
4633         BOOST_UBLAS_INLINE
4634         const_reference operator * () const {
4635             return dereference (iterator_category ());
4636         }
4637         BOOST_UBLAS_INLINE
4638         const_reference operator [] (difference_type n) const {
4639             return *(*this + n);
4640         }
4641
4642         // Index
4643         BOOST_UBLAS_INLINE
4644         size_type index () const {
4645             return it2_.index2 ();
4646         }
4647
4648         // Assignment
4649         BOOST_UBLAS_INLINE
4650         const_iterator &operator = (const const_iterator &it) {
4651             container_const_reference<self_type>::assign (&it ());
4652             it2_ = it.it2_;
4653 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4654             e1_begin_ = it.e1_begin_;
4655             e1_end_ = it.e1_end_;
4656 #endif
4657             return *this;
4658         }
4659
4660         // Comparison
4661         BOOST_UBLAS_INLINE
4662         bool operator == (const const_iterator &it) const {
4663             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
4664             return it2_ == it.it2_;
4665         }
4666         BOOST_UBLAS_INLINE
4667         bool operator < (const const_iterator &it) const {
4668             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
4669             return it2_ < it.it2_;
4670         }
4671
4672     private:
4673         const_subiterator2_type it2_;
4674 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4675         // Mutable due to assignment
4676         /* const */ const_subiterator1_type e1_begin_;
4677         /* const */ const_subiterator1_type e1_end_;
4678 #endif
4679     };
4680 #endif
4681
4682     BOOST_UBLAS_INLINE
4683     const_iterator begin () const {
4684         return find (0);
4685     }
4686     BOOST_UBLAS_INLINE
4687     const_iterator cbegin () const {
4688         return begin ();
4689     }
4690     BOOST_UBLAS_INLINE
4691     const_iterator end () const {
4692         return find (size ());
4693     }
4694     BOOST_UBLAS_INLINE
4695     const_iterator cend () const {
4696         return end ();
4697     }
4698
4699     // Reverse iterator
4700     typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
4701
4702     BOOST_UBLAS_INLINE
4703     const_reverse_iterator rbegin () const {
4704         return const_reverse_iterator (end ());
4705     }
4706     BOOST_UBLAS_INLINE
4707     const_reverse_iterator crbegin () const {
4708         return rbegin ();
4709     }
4710     BOOST_UBLAS_INLINE
4711     const_reverse_iterator rend () const {
4712         return const_reverse_iterator (begin ());
4713     }
4714     BOOST_UBLAS_INLINE
4715     const_reverse_iterator crend () const {
4716         return rend ();
4717     }
4718
4719 private:
4720     expression1_closure_type e1_;
4721     expression2_closure_type e2_;
4722 };
4723
4724 template<class T1, class E1, class T2, class E2>
4725 struct matrix_vector_binary2_traits {
4726     typedef unknown_storage_tag storage_category;
4727     typedef column_major_tag orientation_category;
4728     typedef typename promote_traits<T1, T2>::promote_type promote_type;
4729     typedef matrix_vector_binary2<E1, E2, matrix_vector_prod2<E1, E2, promote_type> > expression_type;
4730 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
4731     typedef expression_type result_type;
4732 #else
4733     typedef typename E2::vector_temporary_type result_type;
4734 #endif
4735 };
4736
4737 template<class E1, class E2>
4738 BOOST_UBLAS_INLINE
4739 typename matrix_vector_binary2_traits<typename E1::value_type, E1,
4740 typename E2::value_type, E2>::result_type
4741 prod (const vector_expression<E1> &e1,
4742       const matrix_expression<E2> &e2,
4743       unknown_storage_tag,
4744       column_major_tag) {
4745     typedef typename matrix_vector_binary2_traits<typename E1::value_type, E1,
4746         typename E2::value_type, E2>::expression_type expression_type;
4747     return expression_type (e1 (), e2 ());
4748 }
4749
4750 // Dispatcher
4751 template<class E1, class E2>
4752 BOOST_UBLAS_INLINE
4753 typename matrix_vector_binary2_traits<typename E1::value_type, E1,
4754 typename E2::value_type, E2>::result_type
4755 prod (const vector_expression<E1> &e1,
4756       const matrix_expression<E2> &e2) {
4757     BOOST_STATIC_ASSERT (E1::complexity == 0);
4758     typedef typename matrix_vector_binary2_traits<typename E1::value_type, E1,
4759         typename E2::value_type, E2>::storage_category storage_category;
4760     typedef typename matrix_vector_binary2_traits<typename E1::value_type, E1,
4761         typename E2::value_type, E2>::orientation_category orientation_category;
4762     return prod (e1, e2, storage_category (), orientation_category ());
4763 }
4764
4765 template<class E1, class E2>
4766 BOOST_UBLAS_INLINE
4767 typename matrix_vector_binary2_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4768 typename type_traits<typename E2::value_type>::precision_type, E2>::result_type
4769 prec_prod (const vector_expression<E1> &e1,
4770            const matrix_expression<E2> &e2,
4771            unknown_storage_tag,
4772            column_major_tag) {
4773     typedef typename matrix_vector_binary2_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4774         typename type_traits<typename E2::value_type>::precision_type, E2>::expression_type expression_type;
4775     return expression_type (e1 (), e2 ());
4776 }
4777
4778 // Dispatcher
4779 template<class E1, class E2>
4780 BOOST_UBLAS_INLINE
4781 typename matrix_vector_binary2_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4782 typename type_traits<typename E2::value_type>::precision_type, E2>::result_type
4783 prec_prod (const vector_expression<E1> &e1,
4784            const matrix_expression<E2> &e2) {
4785     BOOST_STATIC_ASSERT (E1::complexity == 0);
4786     typedef typename matrix_vector_binary2_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4787         typename type_traits<typename E2::value_type>::precision_type, E2>::storage_category storage_category;
4788     typedef typename matrix_vector_binary2_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4789         typename type_traits<typename E2::value_type>::precision_type, E2>::orientation_category orientation_category;
4790     return prec_prod (e1, e2, storage_category (), orientation_category ());
4791 }
4792
4793 template<class V, class E1, class E2>
4794 BOOST_UBLAS_INLINE
4795 V &
4796 prod (const vector_expression<E1> &e1,
4797       const matrix_expression<E2> &e2,
4798       V &v) {
4799     return v.assign (prod (e1, e2));
4800 }
4801
4802 template<class V, class E1, class E2>
4803 BOOST_UBLAS_INLINE
4804 V &
4805 prec_prod (const vector_expression<E1> &e1,
4806            const matrix_expression<E2> &e2,
4807            V &v) {
4808     return v.assign (prec_prod (e1, e2));
4809 }
4810
4811 template<class V, class E1, class E2>
4812 BOOST_UBLAS_INLINE
4813 V
4814 prod (const vector_expression<E1> &e1,
4815       const matrix_expression<E2> &e2) {
4816     return V (prod (e1, e2));
4817 }
4818
4819 template<class V, class E1, class E2>
4820 BOOST_UBLAS_INLINE
4821 V
4822 prec_prod (const vector_expression<E1> &e1,
4823            const matrix_expression<E2> &e2) {
4824     return V (prec_prod (e1, e2));
4825 }
4826
4827 template<class E1, class E2, class F>
4828 class matrix_matrix_binary:
4829     public matrix_expression<matrix_matrix_binary<E1, E2, F> > {
4830
4831 public:
4832     typedef E1 expression1_type;
4833     typedef E2 expression2_type;
4834 private:
4835     typedef F functor_type;
4836 public:
4837     typedef typename E1::const_closure_type expression1_closure_type;
4838     typedef typename E2::const_closure_type expression2_closure_type;
4839 private:
4840     typedef matrix_matrix_binary<E1, E2, F> self_type;
4841 public:
4842 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
4843     using matrix_expression<self_type>::operator ();
4844 #endif
4845     static const unsigned complexity = 1;
4846     typedef typename promote_traits<typename E1::size_type, typename E2::size_type>::promote_type size_type;
4847     typedef typename promote_traits<typename E1::difference_type, typename E2::difference_type>::promote_type difference_type;
4848     typedef typename F::result_type value_type;
4849     typedef value_type const_reference;
4850     typedef const_reference reference;
4851     typedef const self_type const_closure_type;
4852     typedef const_closure_type closure_type;
4853     typedef unknown_orientation_tag orientation_category;
4854     typedef unknown_storage_tag storage_category;
4855
4856     // Construction and destruction
4857     BOOST_UBLAS_INLINE
4858     matrix_matrix_binary (const expression1_type &e1, const expression2_type &e2):
4859       e1_ (e1), e2_ (e2) {}
4860
4861     // Accessors
4862     BOOST_UBLAS_INLINE
4863     size_type size1 () const {
4864         return e1_.size1 ();
4865     }
4866     BOOST_UBLAS_INLINE
4867     size_type size2 () const {
4868         return e2_.size2 ();
4869     }
4870
4871 public:
4872     // Expression accessors
4873     BOOST_UBLAS_INLINE
4874     const expression1_closure_type &expression1 () const {
4875         return e1_;
4876     }
4877     BOOST_UBLAS_INLINE
4878     const expression2_closure_type &expression2 () const {
4879         return e2_;
4880     }
4881
4882 public:
4883     // Element access
4884     BOOST_UBLAS_INLINE
4885     const_reference operator () (size_type i, size_type j) const {
4886         return functor_type::apply (e1_, e2_, i, j);
4887     }
4888
4889     // Element access
4890     BOOST_UBLAS_INLINE
4891     const_reference operator () (size_type i) const {
4892         return functor_type::apply (e1_, e2_, i);
4893     }
4894
4895     // Closure comparison
4896     BOOST_UBLAS_INLINE
4897     bool same_closure (const matrix_matrix_binary &mmb) const {
4898         return (*this).expression1 ().same_closure (mmb.expression1 ()) &&
4899             (*this).expression2 ().same_closure (mmb.expression2 ());
4900     }
4901
4902     // Iterator types
4903 private:
4904     typedef typename E1::const_iterator1 const_iterator11_type;
4905     typedef typename E1::const_iterator2 const_iterator12_type;
4906     typedef typename E2::const_iterator1 const_iterator21_type;
4907     typedef typename E2::const_iterator2 const_iterator22_type;
4908     typedef const value_type *const_pointer;
4909
4910 public:
4911 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
4912     typedef typename iterator_restrict_traits<typename const_iterator11_type::iterator_category,
4913     typename const_iterator22_type::iterator_category>::iterator_category iterator_category;
4914     typedef indexed_const_iterator1<const_closure_type, iterator_category> const_iterator1;
4915     typedef const_iterator1 iterator1;
4916     typedef indexed_const_iterator2<const_closure_type, iterator_category> const_iterator2;
4917     typedef const_iterator2 iterator2;
4918 #else
4919     class const_iterator1;
4920     typedef const_iterator1 iterator1;
4921     class const_iterator2;
4922     typedef const_iterator2 iterator2;
4923 #endif
4924     typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
4925     typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
4926
4927     // Element lookup
4928     BOOST_UBLAS_INLINE
4929     const_iterator1 find1 (int /* rank */, size_type i, size_type j) const {
4930         // FIXME sparse matrix tests fail!
4931         // const_iterator11_type it11 (e1_.find1 (rank, i, 0));
4932         const_iterator11_type it11 (e1_.find1 (0, i, 0));
4933 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
4934         return const_iterator1 (*this, it11.index1 (), j);
4935 #else
4936         // FIXME sparse matrix tests fail!
4937         // const_iterator22_type it22 (e2_.find2 (rank, 0, j));
4938         const_iterator22_type it22 (e2_.find2 (0, 0, j));
4939         return const_iterator1 (*this, it11, it22);
4940 #endif
4941     }
4942     BOOST_UBLAS_INLINE
4943     const_iterator2 find2 (int /* rank */, size_type i, size_type j) const {
4944         // FIXME sparse matrix tests fail!
4945         // const_iterator22_type it22 (e2_.find2 (rank, 0, j));
4946         const_iterator22_type it22 (e2_.find2 (0, 0, j));
4947 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
4948         return const_iterator2 (*this, i, it22.index2 ());
4949 #else
4950         // FIXME sparse matrix tests fail!
4951         // const_iterator11_type it11 (e1_.find1 (rank, i, 0));
4952         const_iterator11_type it11 (e1_.find1 (0, i, 0));
4953         return const_iterator2 (*this, it11, it22);
4954 #endif
4955     }
4956
4957
4958 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
4959     class const_iterator1:
4960         public container_const_reference<matrix_matrix_binary>,
4961         public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
4962         typename E2::const_iterator2::iterator_category>::iterator_category>::template
4963         iterator_base<const_iterator1, value_type>::type {
4964     public:
4965         typedef typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
4966         typename E2::const_iterator2::iterator_category>::iterator_category iterator_category;
4967         typedef typename matrix_matrix_binary::difference_type difference_type;
4968         typedef typename matrix_matrix_binary::value_type value_type;
4969         typedef typename matrix_matrix_binary::const_reference reference;
4970         typedef typename matrix_matrix_binary::const_pointer pointer;
4971
4972         typedef const_iterator2 dual_iterator_type;
4973         typedef const_reverse_iterator2 dual_reverse_iterator_type;
4974
4975         // Construction and destruction
4976 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4977         BOOST_UBLAS_INLINE
4978         const_iterator1 ():
4979           container_const_reference<self_type> (), it1_ (), it2_ (), it2_begin_ (), it2_end_ () {}
4980         BOOST_UBLAS_INLINE
4981         const_iterator1 (const self_type &mmb, const const_iterator11_type &it1, const const_iterator22_type &it2):
4982           container_const_reference<self_type> (mmb), it1_ (it1), it2_ (it2), it2_begin_ (it2.begin ()), it2_end_ (it2.end ()) {}
4983 #else
4984         BOOST_UBLAS_INLINE
4985         const_iterator1 ():
4986           container_const_reference<self_type> (), it1_ (), it2_ () {}
4987         BOOST_UBLAS_INLINE
4988         const_iterator1 (const self_type &mmb, const const_iterator11_type &it1, const const_iterator22_type &it2):
4989           container_const_reference<self_type> (mmb), it1_ (it1), it2_ (it2) {}
4990 #endif
4991
4992     private:
4993         // Random access specialization
4994         BOOST_UBLAS_INLINE
4995         value_type dereference (dense_random_access_iterator_tag) const {
4996             const self_type &mmb = (*this) ();
4997 #ifdef BOOST_UBLAS_USE_INDEXING
4998             return mmb (index1 (), index2 ());
4999 #elif BOOST_UBLAS_USE_ITERATING
5000             difference_type size = BOOST_UBLAS_SAME (mmb.expression1 ().size2 (), mmb.expression2 ().size1 ());
5001 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
5002             return functor_type::apply (size, it1_.begin (), it2_begin_);
5003 #else
5004             return functor_type::apply (size, it1_.begin (), it2_.begin ());
5005 #endif
5006 #else
5007             difference_type size = BOOST_UBLAS_SAME (mmb.expression1 ().size2 (), mmb.expression2 ().size1 ());
5008             if (size >= BOOST_UBLAS_ITERATOR_THRESHOLD)
5009 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
5010                 return functor_type::apply (size, it1_.begin (), it2_begin_);
5011 #else
5012                 return functor_type::apply (size, it1_.begin (), it2_.begin ());
5013 #endif
5014             else
5015                 return mmb (index1 (), index2 ());
5016 #endif
5017         }
5018
5019         // Packed bidirectional specialization
5020         BOOST_UBLAS_INLINE
5021         value_type dereference (packed_random_access_iterator_tag) const {
5022 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
5023             return functor_type::apply (it1_.begin (), it1_.end (),
5024                                         it2_begin_, it2_end_, packed_random_access_iterator_tag ());
5025 #else
5026 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
5027             return functor_type::apply (it1_.begin (), it1_.end (),
5028                                         it2_.begin (), it2_.end (), packed_random_access_iterator_tag ());
5029 #else
5030             return functor_type::apply (boost::numeric::ublas::begin (it1_, iterator1_tag ()),
5031                                         boost::numeric::ublas::end (it1_, iterator1_tag ()),
5032                                         boost::numeric::ublas::begin (it2_, iterator2_tag ()),
5033                                         boost::numeric::ublas::end (it2_, iterator2_tag ()), packed_random_access_iterator_tag ());
5034 #endif
5035 #endif
5036         }
5037
5038         // Sparse bidirectional specialization
5039         BOOST_UBLAS_INLINE
5040         value_type dereference (sparse_bidirectional_iterator_tag) const {
5041 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
5042             return functor_type::apply (it1_.begin (), it1_.end (),
5043                                         it2_begin_, it2_end_, sparse_bidirectional_iterator_tag ());
5044 #else
5045 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
5046             return functor_type::apply (it1_.begin (), it1_.end (),
5047                                         it2_.begin (), it2_.end (), sparse_bidirectional_iterator_tag ());
5048 #else
5049             return functor_type::apply (boost::numeric::ublas::begin (it1_, iterator1_tag ()),
5050                                         boost::numeric::ublas::end (it1_, iterator1_tag ()),
5051                                         boost::numeric::ublas::begin (it2_, iterator2_tag ()),
5052                                         boost::numeric::ublas::end (it2_, iterator2_tag ()), sparse_bidirectional_iterator_tag ());
5053 #endif
5054 #endif
5055         }
5056
5057     public:
5058         // Arithmetic
5059         BOOST_UBLAS_INLINE
5060         const_iterator1 &operator ++ () {
5061             ++ it1_;
5062             return *this;
5063         }
5064         BOOST_UBLAS_INLINE
5065         const_iterator1 &operator -- () {
5066             -- it1_;
5067             return *this;
5068         }
5069         BOOST_UBLAS_INLINE
5070         const_iterator1 &operator += (difference_type n) {
5071             it1_ += n;
5072             return *this;
5073         }
5074         BOOST_UBLAS_INLINE
5075         const_iterator1 &operator -= (difference_type n) {
5076             it1_ -= n;
5077             return *this;
5078         }
5079         BOOST_UBLAS_INLINE
5080         difference_type operator - (const const_iterator1 &it) const {
5081             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
5082             BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
5083             return it1_ - it.it1_;
5084         }
5085
5086         // Dereference
5087         BOOST_UBLAS_INLINE
5088         const_reference operator * () const {
5089             return dereference (iterator_category ());
5090         }
5091         BOOST_UBLAS_INLINE
5092         const_reference operator [] (difference_type n) const {
5093             return *(*this + n);
5094         }
5095
5096 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
5097         BOOST_UBLAS_INLINE
5098 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5099         typename self_type::
5100 #endif
5101         const_iterator2 begin () const {
5102             return (*this) ().find2 (1, index1 (), 0);
5103         }
5104         BOOST_UBLAS_INLINE
5105 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5106         typename self_type::
5107 #endif
5108         const_iterator2 cbegin () const {
5109             return begin ();
5110         }
5111         BOOST_UBLAS_INLINE
5112 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5113         typename self_type::
5114 #endif
5115         const_iterator2 end () const {
5116             return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
5117         }
5118         BOOST_UBLAS_INLINE
5119 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5120         typename self_type::
5121 #endif
5122         const_iterator2 cend () const {
5123             return end ();
5124         }
5125         BOOST_UBLAS_INLINE
5126 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5127         typename self_type::
5128 #endif
5129         const_reverse_iterator2 rbegin () const {
5130             return const_reverse_iterator2 (end ());
5131         }
5132         BOOST_UBLAS_INLINE
5133 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5134         typename self_type::
5135 #endif
5136         const_reverse_iterator2 crbegin () const {
5137             return rbegin ();
5138         }
5139         BOOST_UBLAS_INLINE
5140 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5141         typename self_type::
5142 #endif
5143         const_reverse_iterator2 rend () const {
5144             return const_reverse_iterator2 (begin ());
5145         }
5146         BOOST_UBLAS_INLINE
5147 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5148         typename self_type::
5149 #endif
5150         const_reverse_iterator2 crend () const {
5151             return rend ();
5152         }
5153 #endif
5154
5155         // Indices
5156         BOOST_UBLAS_INLINE
5157         size_type index1 () const {
5158             return it1_.index1 ();
5159         }
5160         BOOST_UBLAS_INLINE
5161         size_type index2 () const {
5162             return it2_.index2 ();
5163         }
5164
5165         // Assignment
5166         BOOST_UBLAS_INLINE
5167         const_iterator1 &operator = (const const_iterator1 &it) {
5168             container_const_reference<self_type>::assign (&it ());
5169             it1_ = it.it1_;
5170             it2_ = it.it2_;
5171 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
5172             it2_begin_ = it.it2_begin_;
5173             it2_end_ = it.it2_end_;
5174 #endif
5175             return *this;
5176         }
5177
5178         // Comparison
5179         BOOST_UBLAS_INLINE
5180         bool operator == (const const_iterator1 &it) const {
5181             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
5182             BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
5183             return it1_ == it.it1_;
5184         }
5185         BOOST_UBLAS_INLINE
5186         bool operator < (const const_iterator1 &it) const {
5187             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
5188             BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
5189             return it1_ < it.it1_;
5190         }
5191
5192     private:
5193         const_iterator11_type it1_;
5194         // Mutable due to assignment
5195         /* const */ const_iterator22_type it2_;
5196 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
5197         /* const */ const_iterator21_type it2_begin_;
5198         /* const */ const_iterator21_type it2_end_;
5199 #endif
5200     };
5201 #endif
5202
5203     BOOST_UBLAS_INLINE
5204     const_iterator1 begin1 () const {
5205         return find1 (0, 0, 0);
5206     }
5207     BOOST_UBLAS_INLINE
5208     const_iterator1 cbegin1 () const {
5209         return begin1 ();
5210     }
5211     BOOST_UBLAS_INLINE
5212     const_iterator1 end1 () const {
5213         return find1 (0, size1 (), 0);
5214     }
5215     BOOST_UBLAS_INLINE
5216     const_iterator1 cend1 () const {
5217         return end1 ();
5218     }
5219
5220 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
5221     class const_iterator2:
5222         public container_const_reference<matrix_matrix_binary>,
5223         public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
5224         typename E2::const_iterator2::iterator_category>::iterator_category>::template
5225         iterator_base<const_iterator2, value_type>::type {
5226     public:
5227         typedef typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
5228         typename E2::const_iterator2::iterator_category>::iterator_category iterator_category;
5229         typedef typename matrix_matrix_binary::difference_type difference_type;
5230         typedef typename matrix_matrix_binary::value_type value_type;
5231         typedef typename matrix_matrix_binary::const_reference reference;
5232         typedef typename matrix_matrix_binary::const_pointer pointer;
5233
5234         typedef const_iterator1 dual_iterator_type;
5235         typedef const_reverse_iterator1 dual_reverse_iterator_type;
5236
5237         // Construction and destruction
5238 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
5239         BOOST_UBLAS_INLINE
5240         const_iterator2 ():
5241           container_const_reference<self_type> (), it1_ (), it2_ (), it1_begin_ (), it1_end_ () {}
5242         BOOST_UBLAS_INLINE
5243         const_iterator2 (const self_type &mmb, const const_iterator11_type &it1, const const_iterator22_type &it2):
5244           container_const_reference<self_type> (mmb), it1_ (it1), it2_ (it2), it1_begin_ (it1.begin ()), it1_end_ (it1.end ()) {}
5245 #else
5246         BOOST_UBLAS_INLINE
5247         const_iterator2 ():
5248           container_const_reference<self_type> (), it1_ (), it2_ () {}
5249         BOOST_UBLAS_INLINE
5250         const_iterator2 (const self_type &mmb, const const_iterator11_type &it1, const const_iterator22_type &it2):
5251           container_const_reference<self_type> (mmb), it1_ (it1), it2_ (it2) {}
5252 #endif
5253
5254     private:
5255         // Random access specialization
5256         BOOST_UBLAS_INLINE
5257         value_type dereference (dense_random_access_iterator_tag) const {
5258             const self_type &mmb = (*this) ();
5259 #ifdef BOOST_UBLAS_USE_INDEXING
5260             return mmb (index1 (), index2 ());
5261 #elif BOOST_UBLAS_USE_ITERATING
5262             difference_type size = BOOST_UBLAS_SAME (mmb.expression1 ().size2 (), mmb.expression2 ().size1 ());
5263 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
5264             return functor_type::apply (size, it1_begin_, it2_.begin ());
5265 #else
5266             return functor_type::apply (size, it1_.begin (), it2_.begin ());
5267 #endif
5268 #else
5269             difference_type size = BOOST_UBLAS_SAME (mmb.expression1 ().size2 (), mmb.expression2 ().size1 ());
5270             if (size >= BOOST_UBLAS_ITERATOR_THRESHOLD)
5271 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
5272                 return functor_type::apply (size, it1_begin_, it2_.begin ());
5273 #else
5274                 return functor_type::apply (size, it1_.begin (), it2_.begin ());
5275 #endif
5276             else
5277                 return mmb (index1 (), index2 ());
5278 #endif
5279         }
5280
5281         // Packed bidirectional specialization
5282         BOOST_UBLAS_INLINE
5283         value_type dereference (packed_random_access_iterator_tag) const {
5284 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
5285             return functor_type::apply (it1_begin_, it1_end_,
5286                                         it2_.begin (), it2_.end (), packed_random_access_iterator_tag ());
5287 #else
5288 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
5289             return functor_type::apply (it1_.begin (), it1_.end (),
5290                                         it2_.begin (), it2_.end (), packed_random_access_iterator_tag ());
5291 #else
5292             return functor_type::apply (boost::numeric::ublas::begin (it1_, iterator1_tag ()),
5293                                         boost::numeric::ublas::end (it1_, iterator1_tag ()),
5294                                         boost::numeric::ublas::begin (it2_, iterator2_tag ()),
5295                                         boost::numeric::ublas::end (it2_, iterator2_tag ()), packed_random_access_iterator_tag ());
5296 #endif
5297 #endif
5298         }
5299
5300         // Sparse bidirectional specialization
5301         BOOST_UBLAS_INLINE
5302         value_type dereference (sparse_bidirectional_iterator_tag) const {
5303 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
5304             return functor_type::apply (it1_begin_, it1_end_,
5305                                         it2_.begin (), it2_.end (), sparse_bidirectional_iterator_tag ());
5306 #else
5307 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
5308             return functor_type::apply (it1_.begin (), it1_.end (),
5309                                         it2_.begin (), it2_.end (), sparse_bidirectional_iterator_tag ());
5310 #else
5311             return functor_type::apply (boost::numeric::ublas::begin (it1_, iterator1_tag ()),
5312                                         boost::numeric::ublas::end (it1_, iterator1_tag ()),
5313                                         boost::numeric::ublas::begin (it2_, iterator2_tag ()),
5314                                         boost::numeric::ublas::end (it2_, iterator2_tag ()), sparse_bidirectional_iterator_tag ());
5315 #endif
5316 #endif
5317         }
5318
5319     public:
5320         // Arithmetic
5321         BOOST_UBLAS_INLINE
5322         const_iterator2 &operator ++ () {
5323             ++ it2_;
5324             return *this;
5325         }
5326         BOOST_UBLAS_INLINE
5327         const_iterator2 &operator -- () {
5328             -- it2_;
5329             return *this;
5330         }
5331         BOOST_UBLAS_INLINE
5332         const_iterator2 &operator += (difference_type n) {
5333             it2_ += n;
5334             return *this;
5335         }
5336         BOOST_UBLAS_INLINE
5337         const_iterator2 &operator -= (difference_type n) {
5338             it2_ -= n;
5339             return *this;
5340         }
5341         BOOST_UBLAS_INLINE
5342         difference_type operator - (const const_iterator2 &it) const {
5343             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
5344             BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
5345             return it2_ - it.it2_;
5346         }
5347
5348         // Dereference
5349         BOOST_UBLAS_INLINE
5350         const_reference operator * () const {
5351             return dereference (iterator_category ());
5352         }
5353         BOOST_UBLAS_INLINE
5354         const_reference operator [] (difference_type n) const {
5355             return *(*this + n);
5356         }
5357
5358 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
5359         BOOST_UBLAS_INLINE
5360 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5361         typename self_type::
5362 #endif
5363         const_iterator1 begin () const {
5364             return (*this) ().find1 (1, 0, index2 ());
5365         }
5366         BOOST_UBLAS_INLINE
5367 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5368         typename self_type::
5369 #endif
5370         const_iterator1 cbegin () const {
5371             return begin ();
5372         }
5373         BOOST_UBLAS_INLINE
5374 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5375         typename self_type::
5376 #endif
5377         const_iterator1 end () const {
5378             return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
5379         }
5380         BOOST_UBLAS_INLINE
5381 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5382         typename self_type::
5383 #endif
5384         const_iterator1 cend () const {
5385             return end ();
5386         }
5387         BOOST_UBLAS_INLINE
5388 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5389         typename self_type::
5390 #endif
5391         const_reverse_iterator1 rbegin () const {
5392             return const_reverse_iterator1 (end ());
5393         }
5394         BOOST_UBLAS_INLINE
5395 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5396         typename self_type::
5397 #endif
5398         const_reverse_iterator1 crbegin () const {
5399             return rbegin ();
5400         }
5401         BOOST_UBLAS_INLINE
5402 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5403         typename self_type::
5404 #endif
5405         const_reverse_iterator1 rend () const {
5406             return const_reverse_iterator1 (begin ());
5407         }
5408         BOOST_UBLAS_INLINE
5409 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5410         typename self_type::
5411 #endif
5412         const_reverse_iterator1 crend () const {
5413             return rend ();
5414         }
5415 #endif
5416
5417         // Indices
5418         BOOST_UBLAS_INLINE
5419         size_type index1 () const {
5420             return it1_.index1 ();
5421         }
5422         BOOST_UBLAS_INLINE
5423         size_type index2 () const {
5424             return it2_.index2 ();
5425         }
5426
5427         // Assignment
5428         BOOST_UBLAS_INLINE
5429         const_iterator2 &operator = (const const_iterator2 &it) {
5430             container_const_reference<self_type>::assign (&it ());
5431             it1_ = it.it1_;
5432             it2_ = it.it2_;
5433 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
5434             it1_begin_ = it.it1_begin_;
5435             it1_end_ = it.it1_end_;
5436 #endif
5437             return *this;
5438         }
5439
5440         // Comparison
5441         BOOST_UBLAS_INLINE
5442         bool operator == (const const_iterator2 &it) const {
5443             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
5444             BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
5445             return it2_ == it.it2_;
5446         }
5447         BOOST_UBLAS_INLINE
5448         bool operator < (const const_iterator2 &it) const {
5449             BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
5450             BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
5451             return it2_ < it.it2_;
5452         }
5453
5454     private:
5455         // Mutable due to assignment
5456         /* const */ const_iterator11_type it1_;
5457         const_iterator22_type it2_;
5458 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
5459         /* const */ const_iterator12_type it1_begin_;
5460         /* const */ const_iterator12_type it1_end_;
5461 #endif
5462     };
5463 #endif
5464
5465     BOOST_UBLAS_INLINE
5466     const_iterator2 begin2 () const {
5467         return find2 (0, 0, 0);
5468     }
5469     BOOST_UBLAS_INLINE
5470     const_iterator2 cbegin2 () const {
5471         return begin2 ();
5472     }
5473     BOOST_UBLAS_INLINE
5474     const_iterator2 end2 () const {
5475         return find2 (0, 0, size2 ());
5476     }
5477     BOOST_UBLAS_INLINE
5478     const_iterator2 cend2 () const {
5479         return end2 ();
5480     }
5481
5482     // Reverse iterators
5483
5484     BOOST_UBLAS_INLINE
5485     const_reverse_iterator1 rbegin1 () const {
5486         return const_reverse_iterator1 (end1 ());
5487     }
5488     BOOST_UBLAS_INLINE
5489     const_reverse_iterator1 crbegin1 () const {
5490         return rbegin1 ();
5491     }
5492     BOOST_UBLAS_INLINE
5493     const_reverse_iterator1 rend1 () const {
5494         return const_reverse_iterator1 (begin1 ());
5495     }
5496     BOOST_UBLAS_INLINE
5497     const_reverse_iterator1 crend1 () const {
5498         return rend1 ();
5499     }
5500
5501     BOOST_UBLAS_INLINE
5502     const_reverse_iterator2 rbegin2 () const {
5503         return const_reverse_iterator2 (end2 ());
5504     }
5505     BOOST_UBLAS_INLINE
5506     const_reverse_iterator2 crbegin2 () const {
5507         return rbegin2 ();
5508     }
5509     BOOST_UBLAS_INLINE
5510     const_reverse_iterator2 rend2 () const {
5511         return const_reverse_iterator2 (begin2 ());
5512     }
5513     BOOST_UBLAS_INLINE
5514     const_reverse_iterator2 crend2 () const {
5515         return rend2 ();
5516     }
5517
5518 private:
5519     expression1_closure_type e1_;
5520     expression2_closure_type e2_;
5521 };
5522
5523 template<class T1, class E1, class T2, class E2>
5524 struct matrix_matrix_binary_traits {
5525     typedef unknown_storage_tag storage_category;
5526     typedef unknown_orientation_tag orientation_category;
5527     typedef typename promote_traits<T1, T2>::promote_type promote_type;
5528     typedef matrix_matrix_binary<E1, E2, matrix_matrix_prod<E1, E2, promote_type> > expression_type;
5529 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
5530     typedef expression_type result_type;
5531 #else
5532     typedef typename E1::matrix_temporary_type result_type;
5533 #endif
5534 };
5535
5536 template<class E1, class E2>
5537 BOOST_UBLAS_INLINE
5538 typename matrix_matrix_binary_traits<typename E1::value_type, E1,
5539 typename E2::value_type, E2>::result_type
5540 prod (const matrix_expression<E1> &e1,
5541       const matrix_expression<E2> &e2,
5542       unknown_storage_tag,
5543       unknown_orientation_tag) {
5544     typedef typename matrix_matrix_binary_traits<typename E1::value_type, E1,
5545         typename E2::value_type, E2>::expression_type expression_type;
5546     return expression_type (e1 (), e2 ());
5547 }
5548
5549 // Dispatcher
5550 template<class E1, class E2>
5551 BOOST_UBLAS_INLINE
5552 typename matrix_matrix_binary_traits<typename E1::value_type, E1,
5553 typename E2::value_type, E2>::result_type
5554 prod (const matrix_expression<E1> &e1,
5555       const matrix_expression<E2> &e2) {
5556     BOOST_STATIC_ASSERT (E1::complexity == 0 && E2::complexity == 0);
5557     typedef typename matrix_matrix_binary_traits<typename E1::value_type, E1,
5558         typename E2::value_type, E2>::storage_category storage_category;
5559     typedef typename matrix_matrix_binary_traits<typename E1::value_type, E1,
5560         typename E2::value_type, E2>::orientation_category orientation_category;
5561     return prod (e1, e2, storage_category (), orientation_category ());
5562 }
5563
5564 template<class E1, class E2>
5565 BOOST_UBLAS_INLINE
5566 typename matrix_matrix_binary_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
5567 typename type_traits<typename E2::value_type>::precision_type, E2>::result_type
5568 prec_prod (const matrix_expression<E1> &e1,
5569            const matrix_expression<E2> &e2,
5570            unknown_storage_tag,
5571            unknown_orientation_tag) {
5572     typedef typename matrix_matrix_binary_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
5573         typename type_traits<typename E2::value_type>::precision_type, E2>::expression_type expression_type;
5574     return expression_type (e1 (), e2 ());
5575 }
5576
5577 // Dispatcher
5578 template<class E1, class E2>
5579 BOOST_UBLAS_INLINE
5580 typename matrix_matrix_binary_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
5581 typename type_traits<typename E2::value_type>::precision_type, E2>::result_type
5582 prec_prod (const matrix_expression<E1> &e1,
5583            const matrix_expression<E2> &e2) {
5584     BOOST_STATIC_ASSERT (E1::complexity == 0 && E2::complexity == 0);
5585     typedef typename matrix_matrix_binary_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
5586         typename type_traits<typename E2::value_type>::precision_type, E2>::storage_category storage_category;
5587     typedef typename matrix_matrix_binary_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
5588         typename type_traits<typename E2::value_type>::precision_type, E2>::orientation_category orientation_category;
5589     return prec_prod (e1, e2, storage_category (), orientation_category ());
5590 }
5591
5592 template<class M, class E1, class E2>
5593 BOOST_UBLAS_INLINE
5594 M &
5595 prod (const matrix_expression<E1> &e1,
5596       const matrix_expression<E2> &e2,
5597       M &m) {
5598     return m.assign (prod (e1, e2));
5599 }
5600
5601 template<class M, class E1, class E2>
5602 BOOST_UBLAS_INLINE
5603 M &
5604 prec_prod (const matrix_expression<E1> &e1,
5605            const matrix_expression<E2> &e2,
5606            M &m) {
5607     return m.assign (prec_prod (e1, e2));
5608 }
5609
5610 template<class M, class E1, class E2>
5611 BOOST_UBLAS_INLINE
5612 M
5613 prod (const matrix_expression<E1> &e1,
5614       const matrix_expression<E2> &e2) {
5615     return M (prod (e1, e2));
5616 }
5617
5618 template<class M, class E1, class E2>
5619 BOOST_UBLAS_INLINE
5620 M
5621 prec_prod (const matrix_expression<E1> &e1,
5622            const matrix_expression<E2> &e2) {
5623     return M (prec_prod (e1, e2));
5624 }
5625
5626 template<class E, class F>
5627 class matrix_scalar_unary:
5628     public scalar_expression<matrix_scalar_unary<E, F> > {
5629 public:
5630     typedef E expression_type;
5631     typedef F functor_type;
5632     typedef typename F::result_type value_type;
5633     typedef typename E::const_closure_type expression_closure_type;
5634
5635     // Construction and destruction
5636     BOOST_UBLAS_INLINE
5637     explicit matrix_scalar_unary (const expression_type &e):
5638       e_ (e) {}
5639
5640 private:
5641     // Expression accessors
5642     BOOST_UBLAS_INLINE
5643     const expression_closure_type &expression () const {
5644         return e_;
5645     }
5646
5647 public:
5648     BOOST_UBLAS_INLINE
5649     operator value_type () const {
5650         return functor_type::apply (e_);
5651     }
5652
5653 private:
5654     expression_closure_type e_;
5655 };
5656
5657 template<class E, class F>
5658 struct matrix_scalar_unary_traits {
5659     typedef matrix_scalar_unary<E, F> expression_type;
5660 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
5661     typedef expression_type result_type;
5662 #else
5663     typedef typename F::result_type result_type;
5664 #endif
5665 };
5666
5667 template<class E>
5668 BOOST_UBLAS_INLINE
5669 typename matrix_scalar_unary_traits<E, matrix_norm_1<E> >::result_type
5670 norm_1 (const matrix_expression<E> &e) {
5671     typedef typename matrix_scalar_unary_traits<E, matrix_norm_1<E> >::expression_type expression_type;
5672     return expression_type (e ());
5673 }
5674
5675 template<class E>
5676 BOOST_UBLAS_INLINE
5677 typename matrix_scalar_unary_traits<E, matrix_norm_frobenius<E> >::result_type
5678 norm_frobenius (const matrix_expression<E> &e) {
5679     typedef typename matrix_scalar_unary_traits<E, matrix_norm_frobenius<E> >::expression_type expression_type;
5680     return expression_type (e ());
5681 }
5682
5683 template<class E>
5684 BOOST_UBLAS_INLINE
5685 typename matrix_scalar_unary_traits<E, matrix_norm_inf<E> >::result_type
5686 norm_inf (const matrix_expression<E> &e) {
5687     typedef typename matrix_scalar_unary_traits<E, matrix_norm_inf<E> >::expression_type expression_type;
5688     return expression_type (e ());
5689 }
5690
5691 }}}
5692
5693 #endif