boost/numeric/ublas/matrix_expression.hpp

   1 //
   2 //  Copyright (c) 2000-2002
   3 //  Joerg Walter, Mathias Koch
   4 //
   5 //  Distributed under the Boost Software License, Version 1.0. (See
   6 //  accompanying file LICENSE_1_0.txt or copy at
   7 //  http://www.boost.org/LICENSE_1_0.txt)
   8 //
   9 //  The authors gratefully acknowledge the support of
  10 //  GeNeSys mbH & Co. KG in producing this work.
  11 //
  12
  13 #ifndef _BOOST_UBLAS_MATRIX_EXPRESSION_
  14 #define _BOOST_UBLAS_MATRIX_EXPRESSION_
  15
  16 #include <boost/numeric/ublas/vector_expression.hpp>
  17
  18 // Expression templates based on ideas of Todd Veldhuizen and Geoffrey Furnish
  19 // Iterators based on ideas of Jeremy Siek
  20 //
  21 // Classes that model the Matrix Expression concept
  22
  23 namespace boost { namespace numeric { namespace ublas {
  24
  25     template<class E>
  26     class matrix_reference:
  27         public matrix_expression<matrix_reference<E> > {
  28
  29         typedef matrix_reference<E> self_type;
  30     public:
  31 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
  32         using matrix_expression<self_type>::operator ();
  33 #endif
  34         typedef typename E::size_type size_type;
  35         typedef typename E::difference_type difference_type;
  36         typedef typename E::value_type value_type;
  37         typedef typename E::const_reference const_reference;
  38         typedef typename boost::mpl::if_<boost::is_const<E>,
  39                                           typename E::const_reference,
  40                                           typename E::reference>::type reference;
  41         typedef E referred_type;
  42         typedef const self_type const_closure_type;
  43         typedef self_type closure_type;
  44         typedef typename E::orientation_category orientation_category;
  45         typedef typename E::storage_category storage_category;
  46
  47         // Construction and destruction
  48         BOOST_UBLAS_INLINE
  49         explicit matrix_reference (referred_type &e):
  50               e_ (e) {}
  51
  52         // Accessors
  53         BOOST_UBLAS_INLINE
  54         size_type size1 () const {
  55             return e_.size1 ();
  56         }
  57         BOOST_UBLAS_INLINE
  58         size_type size2 () const {
  59             return e_.size2 ();
  60         }
  61
  62     public:
  63         // Expression accessors - const correct
  64         BOOST_UBLAS_INLINE
  65         const referred_type &expression () const {
  66             return e_;
  67         }
  68         BOOST_UBLAS_INLINE
  69         referred_type &expression () {
  70             return e_;
  71         }
  72
  73     public:
  74         // Element access
  75 #ifndef BOOST_UBLAS_REFERENCE_CONST_MEMBER
  76         BOOST_UBLAS_INLINE
  77         const_reference operator () (size_type i, size_type j) const {
  78             return expression () (i, j);
  79         }
  80         BOOST_UBLAS_INLINE
  81         reference operator () (size_type i, size_type j) {
  82             return expression () (i, j);
  83         }
  84 #else
  85         BOOST_UBLAS_INLINE
  86         reference operator () (size_type i, size_type j) const {
  87             return expression () (i, j);
  88         }
  89 #endif
  90
  91         // Assignment
  92         BOOST_UBLAS_INLINE
  93         matrix_reference &operator = (const matrix_reference &m) {
  94             expression ().operator = (m);
  95             return *this;
  96         }
  97         template<class AE>
  98         BOOST_UBLAS_INLINE
  99         matrix_reference &operator = (const matrix_expression<AE> &ae) {
 100             expression ().operator = (ae);
 101             return *this;
 102         }
 103         template<class AE>
 104         BOOST_UBLAS_INLINE
 105         matrix_reference &assign (const matrix_expression<AE> &ae) {
 106             expression ().assign (ae);
 107             return *this;
 108         }
 109         template<class AE>
 110         BOOST_UBLAS_INLINE
 111         matrix_reference &operator += (const matrix_expression<AE> &ae) {
 112             expression ().operator += (ae);
 113             return *this;
 114         }
 115         template<class AE>
 116         BOOST_UBLAS_INLINE
 117         matrix_reference &plus_assign (const matrix_expression<AE> &ae) {
 118             expression ().plus_assign (ae);
 119             return *this;
 120         }
 121         template<class AE>
 122         BOOST_UBLAS_INLINE
 123         matrix_reference &operator -= (const matrix_expression<AE> &ae) {
 124             expression ().operator -= (ae);
 125             return *this;
 126         }
 127         template<class AE>
 128         BOOST_UBLAS_INLINE
 129         matrix_reference &minus_assign (const matrix_expression<AE> &ae) {
 130             expression ().minus_assign (ae);
 131             return *this;
 132         }
 133         template<class AT>
 134         BOOST_UBLAS_INLINE
 135         matrix_reference &operator *= (const AT &at) {
 136             expression ().operator *= (at);
 137             return *this;
 138         }
 139         template<class AT>
 140         BOOST_UBLAS_INLINE
 141         matrix_reference &operator /= (const AT &at) {
 142             expression ().operator /= (at);
 143             return *this;
 144         }
 145
 146          // Swapping
 147         BOOST_UBLAS_INLINE
 148         void swap (matrix_reference &m) {
 149             expression ().swap (m.expression ());
 150         }
 151
 152         // Closure comparison
 153         BOOST_UBLAS_INLINE
 154         bool same_closure (const matrix_reference &mr) const {
 155             return &(*this).e_ == &mr.e_;
 156         }
 157
 158         // Iterator types
 159         typedef typename E::const_iterator1 const_iterator1;
 160         typedef typename boost::mpl::if_<boost::is_const<E>,
 161                                           typename E::const_iterator1,
 162                                           typename E::iterator1>::type iterator1;
 163         typedef typename E::const_iterator2 const_iterator2;
 164         typedef typename boost::mpl::if_<boost::is_const<E>,
 165                                           typename E::const_iterator2,
 166                                           typename E::iterator2>::type iterator2;
 167
 168         // Element lookup
 169         BOOST_UBLAS_INLINE
 170         const_iterator1 find1 (int rank, size_type i, size_type j) const {
 171             return expression ().find1 (rank, i, j);
 172         }
 173         BOOST_UBLAS_INLINE
 174         iterator1 find1 (int rank, size_type i, size_type j) {
 175             return expression ().find1 (rank, i, j);
 176         }
 177         BOOST_UBLAS_INLINE
 178         const_iterator2 find2 (int rank, size_type i, size_type j) const {
 179             return expression ().find2 (rank, i, j);
 180         }
 181         BOOST_UBLAS_INLINE
 182         iterator2 find2 (int rank, size_type i, size_type j) {
 183             return expression ().find2 (rank, i, j);
 184         }
 185
 186         // Iterators are the iterators of the referenced expression.
 187
 188         BOOST_UBLAS_INLINE
 189         const_iterator1 begin1 () const {
 190             return expression ().begin1 ();
 191         }
 192         BOOST_UBLAS_INLINE
 193         const_iterator1 cbegin1 () const {
 194             return begin1 ();
 195         }
 196         BOOST_UBLAS_INLINE
 197         const_iterator1 end1 () const {
 198             return expression ().end1 ();
 199         }
 200         BOOST_UBLAS_INLINE
 201         const_iterator1 cend1 () const {
 202             return end1 ();
 203         }
 204
 205         BOOST_UBLAS_INLINE
 206         iterator1 begin1 () {
 207             return expression ().begin1 ();
 208         }
 209         BOOST_UBLAS_INLINE
 210         iterator1 end1 () {
 211             return expression ().end1 ();
 212         }
 213
 214         BOOST_UBLAS_INLINE
 215         const_iterator2 begin2 () const {
 216             return expression ().begin2 ();
 217         }
 218         BOOST_UBLAS_INLINE
 219         const_iterator2 cbegin2 () const {
 220             return begin2 ();
 221         }
 222         BOOST_UBLAS_INLINE
 223         const_iterator2 end2 () const {
 224             return expression ().end2 ();
 225         }
 226         BOOST_UBLAS_INLINE
 227         const_iterator2 cend2 () const {
 228             return end2 ();
 229         }
 230
 231         BOOST_UBLAS_INLINE
 232         iterator2 begin2 () {
 233             return expression ().begin2 ();
 234         }
 235         BOOST_UBLAS_INLINE
 236         iterator2 end2 () {
 237             return expression ().end2 ();
 238         }
 239
 240         // Reverse iterators
 241         typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
 242         typedef reverse_iterator_base1<iterator1> reverse_iterator1;
 243
 244         BOOST_UBLAS_INLINE
 245         const_reverse_iterator1 rbegin1 () const {
 246             return const_reverse_iterator1 (end1 ());
 247         }
 248         BOOST_UBLAS_INLINE
 249         const_reverse_iterator1 crbegin1 () const {
 250             return rbegin1 ();
 251         }
 252         BOOST_UBLAS_INLINE
 253         const_reverse_iterator1 rend1 () const {
 254             return const_reverse_iterator1 (begin1 ());
 255         }
 256         BOOST_UBLAS_INLINE
 257         const_reverse_iterator1 crend1 () const {
 258             return rend1 ();
 259         }
 260
 261         BOOST_UBLAS_INLINE
 262         reverse_iterator1 rbegin1 () {
 263             return reverse_iterator1 (end1 ());
 264         }
 265         BOOST_UBLAS_INLINE
 266         reverse_iterator1 rend1 () {
 267             return reverse_iterator1 (begin1 ());
 268         }
 269
 270         typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
 271         typedef reverse_iterator_base2<iterator2> reverse_iterator2;
 272
 273         BOOST_UBLAS_INLINE
 274         const_reverse_iterator2 rbegin2 () const {
 275             return const_reverse_iterator2 (end2 ());
 276         }
 277         BOOST_UBLAS_INLINE
 278         const_reverse_iterator2 crbegin2 () const {
 279             return rbegin2 ();
 280         }
 281         BOOST_UBLAS_INLINE
 282         const_reverse_iterator2 rend2 () const {
 283             return const_reverse_iterator2 (begin2 ());
 284         }
 285         BOOST_UBLAS_INLINE
 286         const_reverse_iterator2 crend2 () const {
 287             return rend2 ();
 288         }
 289
 290         BOOST_UBLAS_INLINE
 291         reverse_iterator2 rbegin2 () {
 292             return reverse_iterator2 (end2 ());
 293         }
 294         BOOST_UBLAS_INLINE
 295         reverse_iterator2 rend2 () {
 296             return reverse_iterator2 (begin2 ());
 297         }
 298
 299     private:
 300         referred_type &e_;
 301     };
 302
 303
 304     template<class E1, class E2, class F>
 305     class vector_matrix_binary:
 306         public matrix_expression<vector_matrix_binary<E1, E2, F> > {
 307
 308         typedef E1 expression1_type;
 309         typedef E2 expression2_type;
 310     public:
 311         typedef typename E1::const_closure_type expression1_closure_type;
 312         typedef typename E2::const_closure_type expression2_closure_type;
 313     private:
 314         typedef vector_matrix_binary<E1, E2, F> self_type;
 315     public:
 316 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
 317         using matrix_expression<self_type>::operator ();
 318 #endif
 319         typedef F functor_type;
 320         typedef typename promote_traits<typename E1::size_type, typename E2::size_type>::promote_type size_type;
 321         typedef typename promote_traits<typename E1::difference_type, typename E2::difference_type>::promote_type difference_type;
 322         typedef typename F::result_type value_type;
 323         typedef value_type const_reference;
 324         typedef const_reference reference;
 325         typedef const self_type const_closure_type;
 326         typedef const_closure_type closure_type;
 327         typedef unknown_orientation_tag orientation_category;
 328         typedef unknown_storage_tag storage_category;
 329
 330         // Construction and destruction
 331         BOOST_UBLAS_INLINE
 332         vector_matrix_binary (const expression1_type &e1, const expression2_type &e2):
 333             e1_ (e1), e2_ (e2) {}
 334
 335         // Accessors
 336         BOOST_UBLAS_INLINE
 337         size_type size1 () const {
 338             return e1_.size ();
 339         }
 340         BOOST_UBLAS_INLINE
 341         size_type size2 () const {
 342             return e2_.size ();
 343         }
 344
 345     public:
 346         // Expression accessors
 347         BOOST_UBLAS_INLINE
 348         const expression1_closure_type &expression1 () const {
 349             return e1_;
 350         }
 351         BOOST_UBLAS_INLINE
 352         const expression2_closure_type &expression2 () const {
 353             return e2_;
 354         }
 355
 356     public:
 357         // Element access
 358         BOOST_UBLAS_INLINE
 359         const_reference operator () (size_type i, size_type j) const {
 360             return functor_type::apply (e1_ (i), e2_ (j));
 361         }
 362
 363         // Closure comparison
 364         BOOST_UBLAS_INLINE
 365         bool same_closure (const vector_matrix_binary &vmb) const {
 366             return (*this).expression1 ().same_closure (vmb.expression1 ()) &&
 367                    (*this).expression2 ().same_closure (vmb.expression2 ());
 368         }
 369
 370         // Iterator types
 371     private:
 372         typedef typename E1::const_iterator const_subiterator1_type;
 373         typedef typename E2::const_iterator const_subiterator2_type;
 374         typedef const value_type *const_pointer;
 375
 376     public:
 377 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
 378         typedef typename iterator_restrict_traits<typename const_subiterator1_type::iterator_category,
 379                                                   typename const_subiterator2_type::iterator_category>::iterator_category iterator_category;
 380         typedef indexed_const_iterator1<const_closure_type, iterator_category> const_iterator1;
 381         typedef const_iterator1 iterator1;
 382         typedef indexed_const_iterator2<const_closure_type, iterator_category> const_iterator2;
 383         typedef const_iterator2 iterator2;
 384 #else
 385         class const_iterator1;
 386         typedef const_iterator1 iterator1;
 387         class const_iterator2;
 388         typedef const_iterator2 iterator2;
 389 #endif
 390         typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
 391         typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
 392
 393         // Element lookup
 394         BOOST_UBLAS_INLINE
 395         const_iterator1 find1 (int rank, size_type i, size_type j) const {
 396             const_subiterator1_type it1 (e1_.find (i));
 397             const_subiterator1_type it1_end (e1_.find (size1 ()));
 398             const_subiterator2_type it2 (e2_.find (j));
 399             const_subiterator2_type it2_end (e2_.find (size2 ()));
 400             if (it2 == it2_end || (rank == 1 && (it2.index () != j || *it2 == value_type/*zero*/()))) {
 401                 it1 = it1_end;
 402                 it2 = it2_end;
 403             }
 404 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
 405             return const_iterator1 (*this, it1.index (), it2.index ());
 406 #else
 407 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
 408             return const_iterator1 (*this, it1, it2, it2 != it2_end ? *it2 : value_type/*zero*/());
 409 #else
 410             return const_iterator1 (*this, it1, it2);
 411 #endif
 412 #endif
 413         }
 414         BOOST_UBLAS_INLINE
 415         const_iterator2 find2 (int rank, size_type i, size_type j) const {
 416             const_subiterator2_type it2 (e2_.find (j));
 417             const_subiterator2_type it2_end (e2_.find (size2 ()));
 418             const_subiterator1_type it1 (e1_.find (i));
 419             const_subiterator1_type it1_end (e1_.find (size1 ()));
 420             if (it1 == it1_end || (rank == 1 && (it1.index () != i || *it1 == value_type/*zero*/()))) {
 421                 it2 = it2_end;
 422                 it1 = it1_end;
 423             }
 424 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
 425             return const_iterator2 (*this, it1.index (), it2.index ());
 426 #else
 427 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
 428             return const_iterator2 (*this, it1, it2, it1 != it1_end ? *it1 : value_type/*zero*/());
 429 #else
 430             return const_iterator2 (*this, it1, it2);
 431 #endif
 432 #endif
 433         }
 434
 435         // Iterators enhance the iterators of the referenced expressions
 436         // with the binary functor.
 437
 438 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
 439         class const_iterator1:
 440             public container_const_reference<vector_matrix_binary>,
 441             public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
 442                                                                           typename E2::const_iterator::iterator_category>::iterator_category>::template
 443                 iterator_base<const_iterator1, value_type>::type {
 444         public:
 445             typedef typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
 446                                                       typename E2::const_iterator::iterator_category>::iterator_category iterator_category;
 447             typedef typename vector_matrix_binary::difference_type difference_type;
 448             typedef typename vector_matrix_binary::value_type value_type;
 449             typedef typename vector_matrix_binary::const_reference reference;
 450             typedef typename vector_matrix_binary::const_pointer pointer;
 451
 452             typedef const_iterator2 dual_iterator_type;
 453             typedef const_reverse_iterator2 dual_reverse_iterator_type;
 454
 455             // Construction and destruction
 456 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
 457             BOOST_UBLAS_INLINE
 458             const_iterator1 ():
 459                 container_const_reference<self_type> (), it1_ (), it2_ (), t2_ () {}
 460             BOOST_UBLAS_INLINE
 461             const_iterator1 (const self_type &vmb, const const_subiterator1_type &it1, const const_subiterator2_type &it2, value_type t2):
 462                 container_const_reference<self_type> (vmb), it1_ (it1), it2_ (it2), t2_ (t2) {}
 463 #else
 464             BOOST_UBLAS_INLINE
 465             const_iterator1 ():
 466                 container_const_reference<self_type> (), it1_ (), it2_ () {}
 467             BOOST_UBLAS_INLINE
 468             const_iterator1 (const self_type &vmb, const const_subiterator1_type &it1, const const_subiterator2_type &it2):
 469                 container_const_reference<self_type> (vmb), it1_ (it1), it2_ (it2) {}
 470 #endif
 471
 472             // Arithmetic
 473             BOOST_UBLAS_INLINE
 474             const_iterator1 &operator ++ () {
 475                 ++ it1_;
 476                 return *this;
 477             }
 478             BOOST_UBLAS_INLINE
 479             const_iterator1 &operator -- () {
 480                 -- it1_;
 481                 return *this;
 482             }
 483             BOOST_UBLAS_INLINE
 484             const_iterator1 &operator += (difference_type n) {
 485                 it1_ += n;
 486                 return *this;
 487             }
 488             BOOST_UBLAS_INLINE
 489             const_iterator1 &operator -= (difference_type n) {
 490                 it1_ -= n;
 491                 return *this;
 492             }
 493             BOOST_UBLAS_INLINE
 494             difference_type operator - (const const_iterator1 &it) const {
 495                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
 496                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
 497                 return it1_ - it.it1_;
 498             }
 499
 500             // Dereference
 501             BOOST_UBLAS_INLINE
 502             const_reference operator * () const {
 503 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
 504                 return functor_type::apply (*it1_, t2_);
 505 #else
 506                 return functor_type::apply (*it1_, *it2_);
 507 #endif
 508             }
 509             BOOST_UBLAS_INLINE
 510             const_reference operator [] (difference_type n) const {
 511                 return *(*this + n);
 512             }
 513
 514 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
 515             BOOST_UBLAS_INLINE
 516 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 517             typename self_type::
 518 #endif
 519             const_iterator2 begin () const {
 520                 return (*this) ().find2 (1, index1 (), 0);
 521             }
 522             BOOST_UBLAS_INLINE
 523 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 524             typename self_type::
 525 #endif
 526             const_iterator2 cbegin () const {
 527                 return begin ();
 528             }
 529             BOOST_UBLAS_INLINE
 530 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 531             typename self_type::
 532 #endif
 533             const_iterator2 end () const {
 534                 return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
 535             }
 536             BOOST_UBLAS_INLINE
 537 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 538             typename self_type::
 539 #endif
 540             const_iterator2 cend () const {
 541                 return end ();
 542             }
 543             BOOST_UBLAS_INLINE
 544 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 545             typename self_type::
 546 #endif
 547             const_reverse_iterator2 rbegin () const {
 548                 return const_reverse_iterator2 (end ());
 549             }
 550             BOOST_UBLAS_INLINE
 551 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 552             typename self_type::
 553 #endif
 554             const_reverse_iterator2 crbegin () const {
 555                 return rbegin ();
 556             }
 557             BOOST_UBLAS_INLINE
 558 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 559             typename self_type::
 560 #endif
 561             const_reverse_iterator2 rend () const {
 562                 return const_reverse_iterator2 (begin ());
 563             }
 564             BOOST_UBLAS_INLINE
 565 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 566             typename self_type::
 567 #endif
 568             const_reverse_iterator2 crend () const {
 569                 return rend ();
 570             }
 571 #endif
 572
 573             // Indices
 574             BOOST_UBLAS_INLINE
 575             size_type index1 () const {
 576                 return it1_.index ();
 577             }
 578             BOOST_UBLAS_INLINE
 579             size_type  index2 () const {
 580                 return it2_.index ();
 581             }
 582
 583             // Assignment
 584             BOOST_UBLAS_INLINE
 585             const_iterator1 &operator = (const const_iterator1 &it) {
 586                 container_const_reference<self_type>::assign (&it ());
 587                 it1_ = it.it1_;
 588                 it2_ = it.it2_;
 589 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
 590                 t2_ = it.t2_;
 591 #endif
 592                 return *this;
 593             }
 594
 595             // Comparison
 596             BOOST_UBLAS_INLINE
 597             bool operator == (const const_iterator1 &it) const {
 598                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
 599                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
 600                 return it1_ == it.it1_;
 601             }
 602             BOOST_UBLAS_INLINE
 603             bool operator < (const const_iterator1 &it) const {
 604                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
 605                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
 606                 return it1_ < it.it1_;
 607             }
 608
 609         private:
 610 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
 611             const_subiterator1_type it1_;
 612             // Mutable due to assignment
 613             /* const */ const_subiterator2_type it2_;
 614             value_type t2_;
 615 #else
 616             const_subiterator1_type it1_;
 617             const_subiterator2_type it2_;
 618 #endif
 619         };
 620 #endif
 621
 622         BOOST_UBLAS_INLINE
 623         const_iterator1 begin1 () const {
 624             return find1 (0, 0, 0);
 625         }
 626         BOOST_UBLAS_INLINE
 627         const_iterator1 cbegin1 () const {
 628             return begin1 ();
 629         }
 630         BOOST_UBLAS_INLINE
 631         const_iterator1 end1 () const {
 632             return find1 (0, size1 (), 0);
 633         }
 634         BOOST_UBLAS_INLINE
 635         const_iterator1 cend1 () const {
 636             return end1 ();
 637         }
 638
 639 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
 640         class const_iterator2:
 641             public container_const_reference<vector_matrix_binary>,
 642             public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
 643                                                                           typename E2::const_iterator::iterator_category>::iterator_category>::template
 644                 iterator_base<const_iterator2, value_type>::type {
 645         public:
 646             typedef typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
 647                                                       typename E2::const_iterator::iterator_category>::iterator_category iterator_category;
 648             typedef typename vector_matrix_binary::difference_type difference_type;
 649             typedef typename vector_matrix_binary::value_type value_type;
 650             typedef typename vector_matrix_binary::const_reference reference;
 651             typedef typename vector_matrix_binary::const_pointer pointer;
 652
 653             typedef const_iterator1 dual_iterator_type;
 654             typedef const_reverse_iterator1 dual_reverse_iterator_type;
 655
 656             // Construction and destruction
 657 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
 658             BOOST_UBLAS_INLINE
 659             const_iterator2 ():
 660                 container_const_reference<self_type> (), it1_ (), it2_ (), t1_ () {}
 661             BOOST_UBLAS_INLINE
 662             const_iterator2 (const self_type &vmb, const const_subiterator1_type &it1, const const_subiterator2_type &it2, value_type t1):
 663                 container_const_reference<self_type> (vmb), it1_ (it1), it2_ (it2), t1_ (t1) {}
 664 #else
 665             BOOST_UBLAS_INLINE
 666             const_iterator2 ():
 667                 container_const_reference<self_type> (), it1_ (), it2_ () {}
 668             BOOST_UBLAS_INLINE
 669             const_iterator2 (const self_type &vmb, const const_subiterator1_type &it1, const const_subiterator2_type &it2):
 670                 container_const_reference<self_type> (vmb), it1_ (it1), it2_ (it2) {}
 671 #endif
 672
 673             // Arithmetic
 674             BOOST_UBLAS_INLINE
 675             const_iterator2 &operator ++ () {
 676                 ++ it2_;
 677                 return *this;
 678             }
 679             BOOST_UBLAS_INLINE
 680             const_iterator2 &operator -- () {
 681                 -- it2_;
 682                 return *this;
 683             }
 684             BOOST_UBLAS_INLINE
 685             const_iterator2 &operator += (difference_type n) {
 686                 it2_ += n;
 687                 return *this;
 688             }
 689             BOOST_UBLAS_INLINE
 690             const_iterator2 &operator -= (difference_type n) {
 691                 it2_ -= n;
 692                 return *this;
 693             }
 694             BOOST_UBLAS_INLINE
 695             difference_type operator - (const const_iterator2 &it) const {
 696                 BOOST_UBLAS_CHECK ((*this) ().same_closure(it ()), external_logic ());
 697                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
 698                 return it2_ - it.it2_;
 699             }
 700
 701             // Dereference
 702             BOOST_UBLAS_INLINE
 703             const_reference operator * () const {
 704 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
 705                 return functor_type::apply (t1_, *it2_);
 706 #else
 707                 return functor_type::apply (*it1_, *it2_);
 708 #endif
 709             }
 710             BOOST_UBLAS_INLINE
 711             const_reference operator [] (difference_type n) const {
 712                 return *(*this + n);
 713             }
 714
 715 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
 716             BOOST_UBLAS_INLINE
 717 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 718             typename self_type::
 719 #endif
 720             const_iterator1 begin () const {
 721                 return (*this) ().find1 (1, 0, index2 ());
 722             }
 723             BOOST_UBLAS_INLINE
 724 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 725             typename self_type::
 726 #endif
 727             const_iterator1 cbegin () const {
 728                 return begin ();
 729             }
 730             BOOST_UBLAS_INLINE
 731 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 732             typename self_type::
 733 #endif
 734             const_iterator1 end () const {
 735                 return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
 736             }
 737             BOOST_UBLAS_INLINE
 738 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 739             typename self_type::
 740 #endif
 741             const_iterator1 cend () const {
 742                 return end ();
 743             }
 744             BOOST_UBLAS_INLINE
 745 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 746             typename self_type::
 747 #endif
 748             const_reverse_iterator1 rbegin () const {
 749                 return const_reverse_iterator1 (end ());
 750             }
 751             BOOST_UBLAS_INLINE
 752 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 753             typename self_type::
 754 #endif
 755             const_reverse_iterator1 crbegin () const {
 756                 return rbegin ();
 757             }
 758             BOOST_UBLAS_INLINE
 759 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 760             typename self_type::
 761 #endif
 762             const_reverse_iterator1 rend () const {
 763                 return const_reverse_iterator1 (begin ());
 764             }
 765             BOOST_UBLAS_INLINE
 766 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 767             typename self_type::
 768 #endif
 769             const_reverse_iterator1 crend () const {
 770                 return rend ();
 771             }
 772 #endif
 773
 774             // Indices
 775             BOOST_UBLAS_INLINE
 776             size_type index1 () const {
 777                 return it1_.index ();
 778             }
 779             BOOST_UBLAS_INLINE
 780             size_type  index2 () const {
 781                 return it2_.index ();
 782             }
 783
 784             // Assignment
 785             BOOST_UBLAS_INLINE
 786             const_iterator2 &operator = (const const_iterator2 &it) {
 787                 container_const_reference<self_type>::assign (&it ());
 788                 it1_ = it.it1_;
 789                 it2_ = it.it2_;
 790 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
 791                 t1_ = it.t1_;
 792 #endif
 793                 return *this;
 794             }
 795
 796             // Comparison
 797             BOOST_UBLAS_INLINE
 798             bool operator == (const const_iterator2 &it) const {
 799                 BOOST_UBLAS_CHECK ((*this) ().same_closure( it ()), external_logic ());
 800                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
 801                 return it2_ == it.it2_;
 802             }
 803             BOOST_UBLAS_INLINE
 804             bool operator < (const const_iterator2 &it) const {
 805                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
 806                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
 807                 return it2_ < it.it2_;
 808             }
 809
 810         private:
 811 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
 812             // Mutable due to assignment
 813             /* const */ const_subiterator1_type it1_;
 814             const_subiterator2_type it2_;
 815             value_type t1_;
 816 #else
 817             const_subiterator1_type it1_;
 818             const_subiterator2_type it2_;
 819 #endif
 820         };
 821 #endif
 822
 823         BOOST_UBLAS_INLINE
 824         const_iterator2 begin2 () const {
 825             return find2 (0, 0, 0);
 826         }
 827         BOOST_UBLAS_INLINE
 828         const_iterator2 cbegin2 () const {
 829             return begin2 ();
 830         }
 831         BOOST_UBLAS_INLINE
 832         const_iterator2 end2 () const {
 833             return find2 (0, 0, size2 ());
 834         }
 835         BOOST_UBLAS_INLINE
 836         const_iterator2 cend2 () const {
 837             return end2 ();
 838         }
 839
 840         // Reverse iterators
 841
 842         BOOST_UBLAS_INLINE
 843         const_reverse_iterator1 rbegin1 () const {
 844             return const_reverse_iterator1 (end1 ());
 845         }
 846         BOOST_UBLAS_INLINE
 847         const_reverse_iterator1 crbegin1 () const {
 848             return rbegin1 ();
 849         }
 850         BOOST_UBLAS_INLINE
 851         const_reverse_iterator1 rend1 () const {
 852             return const_reverse_iterator1 (begin1 ());
 853         }
 854         BOOST_UBLAS_INLINE
 855         const_reverse_iterator1 crend1 () const {
 856             return rend1 ();
 857         }
 858         BOOST_UBLAS_INLINE
 859         const_reverse_iterator2 rbegin2 () const {
 860             return const_reverse_iterator2 (end2 ());
 861         }
 862         BOOST_UBLAS_INLINE
 863         const_reverse_iterator2 crbegin2 () const {
 864             return rbegin2 ();
 865         }
 866         BOOST_UBLAS_INLINE
 867         const_reverse_iterator2 rend2 () const {
 868             return const_reverse_iterator2 (begin2 ());
 869         }
 870         BOOST_UBLAS_INLINE
 871         const_reverse_iterator2 crend2 () const {
 872             return rend2 ();
 873         }
 874
 875     private:
 876         expression1_closure_type e1_;
 877         expression2_closure_type e2_;
 878     };
 879
 880     template<class E1, class E2, class F>
 881     struct vector_matrix_binary_traits {
 882         typedef vector_matrix_binary<E1, E2, F> expression_type;
 883 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
 884         typedef expression_type result_type;
 885 #else
 886         // ISSUE matrix is arbitary temporary type
 887         typedef matrix<typename F::value_type> result_type;
 888 #endif
 889     };
 890
 891     // (outer_prod (v1, v2)) [i] [j] = v1 [i] * v2 [j]
 892     template<class E1, class E2>
 893     BOOST_UBLAS_INLINE
 894     typename vector_matrix_binary_traits<E1, E2, scalar_multiplies<typename E1::value_type, typename E2::value_type> >::result_type
 895     outer_prod (const vector_expression<E1> &e1,
 896                 const vector_expression<E2> &e2) {
 897         BOOST_STATIC_ASSERT (E1::complexity == 0 && E2::complexity == 0);
 898         typedef typename vector_matrix_binary_traits<E1, E2, scalar_multiplies<typename E1::value_type, typename E2::value_type> >::expression_type expression_type;
 899         return expression_type (e1 (), e2 ());
 900     }
 901
 902     template<class E, class F>
 903     class matrix_unary1:
 904         public matrix_expression<matrix_unary1<E, F> > {
 905
 906         typedef E expression_type;
 907         typedef F functor_type;
 908     public:
 909         typedef typename E::const_closure_type expression_closure_type;
 910     private:
 911         typedef matrix_unary1<E, F> self_type;
 912     public:
 913 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
 914         using matrix_expression<self_type>::operator ();
 915 #endif
 916         typedef typename E::size_type size_type;
 917         typedef typename E::difference_type difference_type;
 918         typedef typename F::result_type value_type;
 919         typedef value_type const_reference;
 920         typedef const_reference reference;
 921         typedef const self_type const_closure_type;
 922         typedef const_closure_type closure_type;
 923         typedef typename E::orientation_category orientation_category;
 924         typedef unknown_storage_tag storage_category;
 925
 926         // Construction and destruction
 927         BOOST_UBLAS_INLINE
 928         explicit matrix_unary1 (const expression_type &e):
 929             e_ (e) {}
 930
 931         // Accessors
 932         BOOST_UBLAS_INLINE
 933         size_type size1 () const {
 934             return e_.size1 ();
 935         }
 936         BOOST_UBLAS_INLINE
 937         size_type size2 () const {
 938             return e_.size2 ();
 939         }
 940
 941     public:
 942         // Expression accessors
 943         BOOST_UBLAS_INLINE
 944         const expression_closure_type &expression () const {
 945             return e_;
 946         }
 947
 948     public:
 949         // Element access
 950         BOOST_UBLAS_INLINE
 951         const_reference operator () (size_type i, size_type j) const {
 952             return functor_type::apply (e_ (i, j));
 953         }
 954
 955         // Closure comparison
 956         BOOST_UBLAS_INLINE
 957         bool same_closure (const matrix_unary1 &mu1) const {
 958             return (*this).expression ().same_closure (mu1.expression ());
 959         }
 960
 961         // Iterator types
 962     private:
 963         typedef typename E::const_iterator1 const_subiterator1_type;
 964         typedef typename E::const_iterator2 const_subiterator2_type;
 965         typedef const value_type *const_pointer;
 966
 967     public:
 968 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
 969         typedef indexed_const_iterator1<const_closure_type, typename const_subiterator1_type::iterator_category> const_iterator1;
 970         typedef const_iterator1 iterator1;
 971         typedef indexed_const_iterator2<const_closure_type, typename const_subiterator2_type::iterator_category> const_iterator2;
 972         typedef const_iterator2 iterator2;
 973 #else
 974         class const_iterator1;
 975         typedef const_iterator1 iterator1;
 976         class const_iterator2;
 977         typedef const_iterator2 iterator2;
 978 #endif
 979         typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
 980         typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
 981
 982         // Element lookup
 983         BOOST_UBLAS_INLINE
 984         const_iterator1 find1 (int rank, size_type i, size_type j) const {
 985             const_subiterator1_type it1 (e_.find1 (rank, i, j));
 986 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
 987             return const_iterator1 (*this, it1.index1 (), it1.index2 ());
 988 #else
 989             return const_iterator1 (*this, it1);
 990 #endif
 991         }
 992         BOOST_UBLAS_INLINE
 993         const_iterator2 find2 (int rank, size_type i, size_type j) const {
 994             const_subiterator2_type it2 (e_.find2 (rank, i, j));
 995 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
 996             return const_iterator2 (*this, it2.index1 (), it2.index2 ());
 997 #else
 998             return const_iterator2 (*this, it2);
 999 #endif
1000         }
1001
1002         // Iterators enhance the iterators of the referenced expression
1003         // with the unary functor.
1004
1005 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1006         class const_iterator1:
1007             public container_const_reference<matrix_unary1>,
1008             public iterator_base_traits<typename E::const_iterator1::iterator_category>::template
1009                 iterator_base<const_iterator1, value_type>::type {
1010         public:
1011             typedef typename E::const_iterator1::iterator_category iterator_category;
1012             typedef typename matrix_unary1::difference_type difference_type;
1013             typedef typename matrix_unary1::value_type value_type;
1014             typedef typename matrix_unary1::const_reference reference;
1015             typedef typename matrix_unary1::const_pointer pointer;
1016
1017             typedef const_iterator2 dual_iterator_type;
1018             typedef const_reverse_iterator2 dual_reverse_iterator_type;
1019
1020             // Construction and destruction
1021             BOOST_UBLAS_INLINE
1022             const_iterator1 ():
1023                 container_const_reference<self_type> (), it_ () {}
1024             BOOST_UBLAS_INLINE
1025             const_iterator1 (const self_type &mu, const const_subiterator1_type &it):
1026                 container_const_reference<self_type> (mu), it_ (it) {}
1027
1028             // Arithmetic
1029             BOOST_UBLAS_INLINE
1030             const_iterator1 &operator ++ () {
1031                 ++ it_;
1032                 return *this;
1033             }
1034             BOOST_UBLAS_INLINE
1035             const_iterator1 &operator -- () {
1036                 -- it_;
1037                 return *this;
1038             }
1039             BOOST_UBLAS_INLINE
1040             const_iterator1 &operator += (difference_type n) {
1041                 it_ += n;
1042                 return *this;
1043             }
1044             BOOST_UBLAS_INLINE
1045             const_iterator1 &operator -= (difference_type n) {
1046                 it_ -= n;
1047                 return *this;
1048             }
1049             BOOST_UBLAS_INLINE
1050             difference_type operator - (const const_iterator1 &it) const {
1051                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1052                 return it_ - it.it_;
1053             }
1054
1055             // Dereference
1056             BOOST_UBLAS_INLINE
1057             const_reference operator * () const {
1058                 return functor_type::apply (*it_);
1059             }
1060             BOOST_UBLAS_INLINE
1061             const_reference operator [] (difference_type n) const {
1062                 return *(*this + n);
1063             }
1064
1065 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1066             BOOST_UBLAS_INLINE
1067 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1068             typename self_type::
1069 #endif
1070             const_iterator2 begin () const {
1071                 return (*this) ().find2 (1, index1 (), 0);
1072             }
1073             BOOST_UBLAS_INLINE
1074 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1075             typename self_type::
1076 #endif
1077             const_iterator2 cbegin () const {
1078                 return begin ();
1079             }
1080             BOOST_UBLAS_INLINE
1081 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1082             typename self_type::
1083 #endif
1084             const_iterator2 end () const {
1085                 return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
1086             }
1087             BOOST_UBLAS_INLINE
1088 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1089             typename self_type::
1090 #endif
1091             const_iterator2 cend () const {
1092                 return end ();
1093             }
1094             BOOST_UBLAS_INLINE
1095 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1096             typename self_type::
1097 #endif
1098             const_reverse_iterator2 rbegin () const {
1099                 return const_reverse_iterator2 (end ());
1100             }
1101             BOOST_UBLAS_INLINE
1102 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1103             typename self_type::
1104 #endif
1105             const_reverse_iterator2 crbegin () const {
1106                 return rbegin ();
1107             }
1108             BOOST_UBLAS_INLINE
1109 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1110             typename self_type::
1111 #endif
1112             const_reverse_iterator2 rend () const {
1113                 return const_reverse_iterator2 (begin ());
1114             }
1115             BOOST_UBLAS_INLINE
1116 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1117             typename self_type::
1118 #endif
1119             const_reverse_iterator2 crend () const {
1120                 return rend ();
1121             }
1122 #endif
1123
1124             // Indices
1125             BOOST_UBLAS_INLINE
1126             size_type index1 () const {
1127                 return it_.index1 ();
1128             }
1129             BOOST_UBLAS_INLINE
1130             size_type index2 () const {
1131                 return it_.index2 ();
1132             }
1133
1134             // Assignment
1135             BOOST_UBLAS_INLINE
1136             const_iterator1 &operator = (const const_iterator1 &it) {
1137                 container_const_reference<self_type>::assign (&it ());
1138                 it_ = it.it_;
1139                 return *this;
1140             }
1141
1142             // Comparison
1143             BOOST_UBLAS_INLINE
1144             bool operator == (const const_iterator1 &it) const {
1145                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1146                 return it_ == it.it_;
1147             }
1148             BOOST_UBLAS_INLINE
1149             bool operator < (const const_iterator1 &it) const {
1150                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1151                 return it_ < it.it_;
1152             }
1153
1154         private:
1155             const_subiterator1_type it_;
1156         };
1157 #endif
1158
1159         BOOST_UBLAS_INLINE
1160         const_iterator1 begin1 () const {
1161             return find1 (0, 0, 0);
1162         }
1163         BOOST_UBLAS_INLINE
1164         const_iterator1 cbegin1 () const {
1165             return begin1 ();
1166         }
1167         BOOST_UBLAS_INLINE
1168         const_iterator1 end1 () const {
1169             return find1 (0, size1 (), 0);
1170         }
1171         BOOST_UBLAS_INLINE
1172         const_iterator1 cend1 () const {
1173             return end1 ();
1174         }
1175
1176 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1177         class const_iterator2:
1178             public container_const_reference<matrix_unary1>,
1179             public iterator_base_traits<typename E::const_iterator2::iterator_category>::template
1180                 iterator_base<const_iterator2, value_type>::type {
1181         public:
1182             typedef typename E::const_iterator2::iterator_category iterator_category;
1183             typedef typename matrix_unary1::difference_type difference_type;
1184             typedef typename matrix_unary1::value_type value_type;
1185             typedef typename matrix_unary1::const_reference reference;
1186             typedef typename matrix_unary1::const_pointer pointer;
1187
1188             typedef const_iterator1 dual_iterator_type;
1189             typedef const_reverse_iterator1 dual_reverse_iterator_type;
1190
1191             // Construction and destruction
1192             BOOST_UBLAS_INLINE
1193             const_iterator2 ():
1194                 container_const_reference<self_type> (), it_ () {}
1195             BOOST_UBLAS_INLINE
1196             const_iterator2 (const self_type &mu, const const_subiterator2_type &it):
1197                 container_const_reference<self_type> (mu), it_ (it) {}
1198
1199             // Arithmetic
1200             BOOST_UBLAS_INLINE
1201             const_iterator2 &operator ++ () {
1202                 ++ it_;
1203                 return *this;
1204             }
1205             BOOST_UBLAS_INLINE
1206             const_iterator2 &operator -- () {
1207                 -- it_;
1208                 return *this;
1209             }
1210             BOOST_UBLAS_INLINE
1211             const_iterator2 &operator += (difference_type n) {
1212                 it_ += n;
1213                 return *this;
1214             }
1215             BOOST_UBLAS_INLINE
1216             const_iterator2 &operator -= (difference_type n) {
1217                 it_ -= n;
1218                 return *this;
1219             }
1220             BOOST_UBLAS_INLINE
1221             difference_type operator - (const const_iterator2 &it) const {
1222                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1223                 return it_ - it.it_;
1224             }
1225
1226             // Dereference
1227             BOOST_UBLAS_INLINE
1228             const_reference operator * () const {
1229                 return functor_type::apply (*it_);
1230             }
1231             BOOST_UBLAS_INLINE
1232             const_reference operator [] (difference_type n) const {
1233                 return *(*this + n);
1234             }
1235
1236 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1237             BOOST_UBLAS_INLINE
1238 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1239             typename self_type::
1240 #endif
1241             const_iterator1 begin () const {
1242                 return (*this) ().find1 (1, 0, index2 ());
1243             }
1244             BOOST_UBLAS_INLINE
1245 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1246             typename self_type::
1247 #endif
1248             const_iterator1 cbegin () const {
1249                 return begin ();
1250             }
1251             BOOST_UBLAS_INLINE
1252 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1253             typename self_type::
1254 #endif
1255             const_iterator1 end () const {
1256                 return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
1257             }
1258             BOOST_UBLAS_INLINE
1259 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1260             typename self_type::
1261 #endif
1262             const_iterator1 cend () const {
1263                 return end ();
1264             }
1265             BOOST_UBLAS_INLINE
1266 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1267             typename self_type::
1268 #endif
1269             const_reverse_iterator1 rbegin () const {
1270                 return const_reverse_iterator1 (end ());
1271             }
1272             BOOST_UBLAS_INLINE
1273 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1274             typename self_type::
1275 #endif
1276             const_reverse_iterator1 crbegin () const {
1277                 return rbegin ();
1278             }
1279             BOOST_UBLAS_INLINE
1280 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1281             typename self_type::
1282 #endif
1283             const_reverse_iterator1 rend () const {
1284                 return const_reverse_iterator1 (begin ());
1285             }
1286             BOOST_UBLAS_INLINE
1287 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1288             typename self_type::
1289 #endif
1290             const_reverse_iterator1 crend () const {
1291                 return rend ();
1292             }
1293 #endif
1294
1295             // Indices
1296             BOOST_UBLAS_INLINE
1297             size_type index1 () const {
1298                 return it_.index1 ();
1299             }
1300             BOOST_UBLAS_INLINE
1301             size_type index2 () const {
1302                 return it_.index2 ();
1303             }
1304
1305             // Assignment
1306             BOOST_UBLAS_INLINE
1307             const_iterator2 &operator = (const const_iterator2 &it) {
1308                 container_const_reference<self_type>::assign (&it ());
1309                 it_ = it.it_;
1310                 return *this;
1311             }
1312
1313             // Comparison
1314             BOOST_UBLAS_INLINE
1315             bool operator == (const const_iterator2 &it) const {
1316                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1317                 return it_ == it.it_;
1318             }
1319             BOOST_UBLAS_INLINE
1320             bool operator < (const const_iterator2 &it) const {
1321                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1322                 return it_ < it.it_;
1323             }
1324
1325         private:
1326             const_subiterator2_type it_;
1327         };
1328 #endif
1329
1330         BOOST_UBLAS_INLINE
1331         const_iterator2 begin2 () const {
1332             return find2 (0, 0, 0);
1333         }
1334         BOOST_UBLAS_INLINE
1335         const_iterator2 cbegin2 () const {
1336             return begin2 ();
1337         }
1338         BOOST_UBLAS_INLINE
1339         const_iterator2 end2 () const {
1340             return find2 (0, 0, size2 ());
1341         }
1342         BOOST_UBLAS_INLINE
1343         const_iterator2 cend2 () const {
1344             return end2 ();
1345         }
1346
1347         // Reverse iterators
1348
1349         BOOST_UBLAS_INLINE
1350         const_reverse_iterator1 rbegin1 () const {
1351             return const_reverse_iterator1 (end1 ());
1352         }
1353         BOOST_UBLAS_INLINE
1354         const_reverse_iterator1 crbegin1 () const {
1355             return rbegin1 ();
1356         }
1357         BOOST_UBLAS_INLINE
1358         const_reverse_iterator1 rend1 () const {
1359             return const_reverse_iterator1 (begin1 ());
1360         }
1361         BOOST_UBLAS_INLINE
1362         const_reverse_iterator1 crend1 () const {
1363             return rend1 ();
1364         }
1365
1366         BOOST_UBLAS_INLINE
1367         const_reverse_iterator2 rbegin2 () const {
1368             return const_reverse_iterator2 (end2 ());
1369         }
1370         BOOST_UBLAS_INLINE
1371         const_reverse_iterator2 crbegin2 () const {
1372             return rbegin2 ();
1373         }
1374         BOOST_UBLAS_INLINE
1375         const_reverse_iterator2 rend2 () const {
1376             return const_reverse_iterator2 (begin2 ());
1377         }
1378         BOOST_UBLAS_INLINE
1379         const_reverse_iterator2 crend2 () const {
1380             return rend2 ();
1381         }
1382
1383     private:
1384         expression_closure_type e_;
1385     };
1386
1387     template<class E, class F>
1388     struct matrix_unary1_traits {
1389         typedef matrix_unary1<E, F> expression_type;
1390 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
1391         typedef expression_type result_type;
1392 #else
1393         typedef typename E::matrix_temporary_type result_type;
1394 #endif
1395     };
1396
1397     // (- m) [i] [j] = - m [i] [j]
1398     template<class E>
1399     BOOST_UBLAS_INLINE
1400     typename matrix_unary1_traits<E, scalar_negate<typename E::value_type> >::result_type
1401     operator - (const matrix_expression<E> &e) {
1402         typedef typename matrix_unary1_traits<E, scalar_negate<typename E::value_type> >::expression_type expression_type;
1403         return expression_type (e ());
1404     }
1405
1406     // (conj m) [i] [j] = conj (m [i] [j])
1407     template<class E>
1408     BOOST_UBLAS_INLINE
1409     typename matrix_unary1_traits<E, scalar_conj<typename E::value_type> >::result_type
1410     conj (const matrix_expression<E> &e) {
1411         typedef typename matrix_unary1_traits<E, scalar_conj<typename E::value_type> >::expression_type expression_type;
1412         return expression_type (e ());
1413     }
1414
1415     // (real m) [i] [j] = real (m [i] [j])
1416     template<class E>
1417     BOOST_UBLAS_INLINE
1418     typename matrix_unary1_traits<E, scalar_real<typename E::value_type> >::result_type
1419     real (const matrix_expression<E> &e) {
1420         typedef typename matrix_unary1_traits<E, scalar_real<typename E::value_type> >::expression_type expression_type;
1421         return expression_type (e ());
1422     }
1423
1424     // (imag m) [i] [j] = imag (m [i] [j])
1425     template<class E>
1426     BOOST_UBLAS_INLINE
1427     typename matrix_unary1_traits<E, scalar_imag<typename E::value_type> >::result_type
1428     imag (const matrix_expression<E> &e) {
1429         typedef typename matrix_unary1_traits<E, scalar_imag<typename E::value_type> >::expression_type expression_type;
1430         return expression_type (e ());
1431     }
1432
1433     template<class E, class F>
1434     class matrix_unary2:
1435         public matrix_expression<matrix_unary2<E, F> > {
1436
1437         typedef typename boost::mpl::if_<boost::is_same<F, scalar_identity<typename E::value_type> >,
1438                                           E,
1439                                           const E>::type expression_type;
1440         typedef F functor_type;
1441     public:
1442         typedef typename boost::mpl::if_<boost::is_const<expression_type>,
1443                                           typename E::const_closure_type,
1444                                           typename E::closure_type>::type expression_closure_type;
1445     private:
1446         typedef matrix_unary2<E, F> self_type;
1447     public:
1448 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
1449         using matrix_expression<self_type>::operator ();
1450 #endif
1451         typedef typename E::size_type size_type;
1452         typedef typename E::difference_type difference_type;
1453         typedef typename F::result_type value_type;
1454         typedef value_type const_reference;
1455         typedef typename boost::mpl::if_<boost::is_same<F, scalar_identity<value_type> >,
1456                                           typename E::reference,
1457                                           value_type>::type reference;
1458
1459         typedef const self_type const_closure_type;
1460         typedef self_type closure_type;
1461         typedef typename boost::mpl::if_<boost::is_same<typename E::orientation_category,
1462                                                          row_major_tag>,
1463                                           column_major_tag,
1464                 typename boost::mpl::if_<boost::is_same<typename E::orientation_category,
1465                                                          column_major_tag>,
1466                                           row_major_tag,
1467                                           typename E::orientation_category>::type>::type orientation_category;
1468         typedef typename E::storage_category storage_category;
1469
1470         // Construction and destruction
1471         BOOST_UBLAS_INLINE
1472         // matrix_unary2 may be used as mutable expression -
1473         // this is the only non const expression constructor
1474         explicit matrix_unary2 (expression_type &e):
1475             e_ (e) {}
1476
1477         // Accessors
1478         BOOST_UBLAS_INLINE
1479         size_type size1 () const {
1480             return e_.size2 ();
1481         }
1482         BOOST_UBLAS_INLINE
1483         size_type size2 () const {
1484             return e_.size1 ();
1485         }
1486
1487     public:
1488         // Expression accessors
1489         BOOST_UBLAS_INLINE
1490         const expression_closure_type &expression () const {
1491             return e_;
1492         }
1493
1494     public:
1495         // Element access
1496         BOOST_UBLAS_INLINE
1497         const_reference operator () (size_type i, size_type j) const {
1498             return functor_type::apply (e_ (j, i));
1499         }
1500         BOOST_UBLAS_INLINE
1501         reference operator () (size_type i, size_type j) {
1502             BOOST_STATIC_ASSERT ((boost::is_same<functor_type, scalar_identity<value_type > >::value));
1503             return e_ (j, i);
1504         }
1505
1506         // Closure comparison
1507         BOOST_UBLAS_INLINE
1508         bool same_closure (const matrix_unary2 &mu2) const {
1509             return (*this).expression ().same_closure (mu2.expression ());
1510         }
1511
1512         // Iterator types
1513     private:
1514         typedef typename E::const_iterator1 const_subiterator2_type;
1515         typedef typename E::const_iterator2 const_subiterator1_type;
1516         typedef const value_type *const_pointer;
1517
1518     public:
1519 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1520         typedef indexed_const_iterator1<const_closure_type, typename const_subiterator1_type::iterator_category> const_iterator1;
1521         typedef const_iterator1 iterator1;
1522         typedef indexed_const_iterator2<const_closure_type, typename const_subiterator2_type::iterator_category> const_iterator2;
1523         typedef const_iterator2 iterator2;
1524 #else
1525         class const_iterator1;
1526         typedef const_iterator1 iterator1;
1527         class const_iterator2;
1528         typedef const_iterator2 iterator2;
1529 #endif
1530         typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
1531         typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
1532
1533         // Element lookup
1534         BOOST_UBLAS_INLINE
1535         const_iterator1 find1 (int rank, size_type i, size_type j) const {
1536             const_subiterator1_type it1 (e_.find2 (rank, j, i));
1537 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1538             return const_iterator1 (*this, it1.index2 (), it1.index1 ());
1539 #else
1540             return const_iterator1 (*this, it1);
1541 #endif
1542         }
1543         BOOST_UBLAS_INLINE
1544         const_iterator2 find2 (int rank, size_type i, size_type j) const {
1545             const_subiterator2_type it2 (e_.find1 (rank, j, i));
1546 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1547             return const_iterator2 (*this, it2.index2 (), it2.index1 ());
1548 #else
1549             return const_iterator2 (*this, it2);
1550 #endif
1551         }
1552
1553         // Iterators enhance the iterators of the referenced expression
1554         // with the unary functor.
1555
1556 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1557         class const_iterator1:
1558             public container_const_reference<matrix_unary2>,
1559             public iterator_base_traits<typename E::const_iterator2::iterator_category>::template
1560                 iterator_base<const_iterator1, value_type>::type {
1561         public:
1562             typedef typename E::const_iterator2::iterator_category iterator_category;
1563             typedef typename matrix_unary2::difference_type difference_type;
1564             typedef typename matrix_unary2::value_type value_type;
1565             typedef typename matrix_unary2::const_reference reference;
1566             typedef typename matrix_unary2::const_pointer pointer;
1567
1568             typedef const_iterator2 dual_iterator_type;
1569             typedef const_reverse_iterator2 dual_reverse_iterator_type;
1570
1571             // Construction and destruction
1572             BOOST_UBLAS_INLINE
1573             const_iterator1 ():
1574                 container_const_reference<self_type> (), it_ () {}
1575             BOOST_UBLAS_INLINE
1576             const_iterator1 (const self_type &mu, const const_subiterator1_type &it):
1577                 container_const_reference<self_type> (mu), it_ (it) {}
1578
1579             // Arithmetic
1580             BOOST_UBLAS_INLINE
1581             const_iterator1 &operator ++ () {
1582                 ++ it_;
1583                 return *this;
1584             }
1585             BOOST_UBLAS_INLINE
1586             const_iterator1 &operator -- () {
1587                 -- it_;
1588                 return *this;
1589             }
1590             BOOST_UBLAS_INLINE
1591             const_iterator1 &operator += (difference_type n) {
1592                 it_ += n;
1593                 return *this;
1594             }
1595             BOOST_UBLAS_INLINE
1596             const_iterator1 &operator -= (difference_type n) {
1597                 it_ -= n;
1598                 return *this;
1599             }
1600             BOOST_UBLAS_INLINE
1601             difference_type operator - (const const_iterator1 &it) const {
1602                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1603                 return it_ - it.it_;
1604             }
1605
1606             // Dereference
1607             BOOST_UBLAS_INLINE
1608             const_reference operator * () const {
1609                 return functor_type::apply (*it_);
1610             }
1611             BOOST_UBLAS_INLINE
1612             const_reference operator [] (difference_type n) const {
1613                 return *(*this + n);
1614             }
1615
1616 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1617             BOOST_UBLAS_INLINE
1618 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1619             typename self_type::
1620 #endif
1621             const_iterator2 begin () const {
1622                 return (*this) ().find2 (1, index1 (), 0);
1623             }
1624             BOOST_UBLAS_INLINE
1625 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1626             typename self_type::
1627 #endif
1628             const_iterator2 cbegin () const {
1629                 return begin ();
1630             }
1631             BOOST_UBLAS_INLINE
1632 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1633             typename self_type::
1634 #endif
1635             const_iterator2 end () const {
1636                 return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
1637             }
1638             BOOST_UBLAS_INLINE
1639 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1640             typename self_type::
1641 #endif
1642             const_iterator2 cend () const {
1643                 return end ();
1644             }
1645             BOOST_UBLAS_INLINE
1646 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1647             typename self_type::
1648 #endif
1649             const_reverse_iterator2 rbegin () const {
1650                 return const_reverse_iterator2 (end ());
1651             }
1652             BOOST_UBLAS_INLINE
1653 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1654             typename self_type::
1655 #endif
1656             const_reverse_iterator2 crbegin () const {
1657                 return rbegin ();
1658             }
1659             BOOST_UBLAS_INLINE
1660 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1661             typename self_type::
1662 #endif
1663             const_reverse_iterator2 rend () const {
1664                 return const_reverse_iterator2 (begin ());
1665             }
1666             BOOST_UBLAS_INLINE
1667 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1668             typename self_type::
1669 #endif
1670             const_reverse_iterator2 crend () const {
1671                 return rend ();
1672             }
1673 #endif
1674
1675             // Indices
1676             BOOST_UBLAS_INLINE
1677             size_type index1 () const {
1678                 return it_.index2 ();
1679             }
1680             BOOST_UBLAS_INLINE
1681             size_type index2 () const {
1682                 return it_.index1 ();
1683             }
1684
1685             // Assignment
1686             BOOST_UBLAS_INLINE
1687             const_iterator1 &operator = (const const_iterator1 &it) {
1688                 container_const_reference<self_type>::assign (&it ());
1689                 it_ = it.it_;
1690                 return *this;
1691             }
1692
1693             // Comparison
1694             BOOST_UBLAS_INLINE
1695             bool operator == (const const_iterator1 &it) const {
1696                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1697                 return it_ == it.it_;
1698             }
1699             BOOST_UBLAS_INLINE
1700             bool operator < (const const_iterator1 &it) const {
1701                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1702                 return it_ < it.it_;
1703             }
1704
1705         private:
1706             const_subiterator1_type it_;
1707         };
1708 #endif
1709
1710         BOOST_UBLAS_INLINE
1711         const_iterator1 begin1 () const {
1712             return find1 (0, 0, 0);
1713         }
1714         BOOST_UBLAS_INLINE
1715         const_iterator1 cbegin1 () const {
1716             return begin1 ();
1717         }
1718         BOOST_UBLAS_INLINE
1719         const_iterator1 end1 () const {
1720             return find1 (0, size1 (), 0);
1721         }
1722         BOOST_UBLAS_INLINE
1723         const_iterator1 cend1 () const {
1724             return end1 ();
1725         }
1726
1727 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1728         class const_iterator2:
1729             public container_const_reference<matrix_unary2>,
1730             public iterator_base_traits<typename E::const_iterator1::iterator_category>::template
1731                 iterator_base<const_iterator2, value_type>::type {
1732         public:
1733             typedef typename E::const_iterator1::iterator_category iterator_category;
1734             typedef typename matrix_unary2::difference_type difference_type;
1735             typedef typename matrix_unary2::value_type value_type;
1736             typedef typename matrix_unary2::const_reference reference;
1737             typedef typename matrix_unary2::const_pointer pointer;
1738
1739             typedef const_iterator1 dual_iterator_type;
1740             typedef const_reverse_iterator1 dual_reverse_iterator_type;
1741
1742             // Construction and destruction
1743             BOOST_UBLAS_INLINE
1744             const_iterator2 ():
1745                 container_const_reference<self_type> (), it_ () {}
1746             BOOST_UBLAS_INLINE
1747             const_iterator2 (const self_type &mu, const const_subiterator2_type &it):
1748                 container_const_reference<self_type> (mu), it_ (it) {}
1749
1750             // Arithmetic
1751             BOOST_UBLAS_INLINE
1752             const_iterator2 &operator ++ () {
1753                 ++ it_;
1754                 return *this;
1755             }
1756             BOOST_UBLAS_INLINE
1757             const_iterator2 &operator -- () {
1758                 -- it_;
1759                 return *this;
1760             }
1761             BOOST_UBLAS_INLINE
1762             const_iterator2 &operator += (difference_type n) {
1763                 it_ += n;
1764                 return *this;
1765             }
1766             BOOST_UBLAS_INLINE
1767             const_iterator2 &operator -= (difference_type n) {
1768                 it_ -= n;
1769                 return *this;
1770             }
1771             BOOST_UBLAS_INLINE
1772             difference_type operator - (const const_iterator2 &it) const {
1773                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1774                 return it_ - it.it_;
1775             }
1776
1777             // Dereference
1778             BOOST_UBLAS_INLINE
1779             const_reference operator * () const {
1780                 return functor_type::apply (*it_);
1781             }
1782             BOOST_UBLAS_INLINE
1783             const_reference operator [] (difference_type n) const {
1784                 return *(*this + n);
1785             }
1786
1787 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1788             BOOST_UBLAS_INLINE
1789 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1790             typename self_type::
1791 #endif
1792             const_iterator1 begin () const {
1793                 return (*this) ().find1 (1, 0, index2 ());
1794             }
1795             BOOST_UBLAS_INLINE
1796 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1797             typename self_type::
1798 #endif
1799             const_iterator1 cbegin () const {
1800                 return begin ();
1801             }
1802             BOOST_UBLAS_INLINE
1803 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1804             typename self_type::
1805 #endif
1806             const_iterator1 end () const {
1807                 return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
1808             }
1809             BOOST_UBLAS_INLINE
1810 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1811             typename self_type::
1812 #endif
1813             const_iterator1 cend () const {
1814                 return end ();
1815             }
1816             BOOST_UBLAS_INLINE
1817 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1818             typename self_type::
1819 #endif
1820             const_reverse_iterator1 rbegin () const {
1821                 return const_reverse_iterator1 (end ());
1822             }
1823             BOOST_UBLAS_INLINE
1824 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1825             typename self_type::
1826 #endif
1827             const_reverse_iterator1 crbegin () const {
1828                 return rbegin ();
1829             }
1830             BOOST_UBLAS_INLINE
1831 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1832             typename self_type::
1833 #endif
1834             const_reverse_iterator1 rend () const {
1835                 return const_reverse_iterator1 (begin ());
1836             }
1837             BOOST_UBLAS_INLINE
1838 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1839             typename self_type::
1840 #endif
1841             const_reverse_iterator1 crend () const {
1842                 return rend ();
1843             }
1844 #endif
1845
1846             // Indices
1847             BOOST_UBLAS_INLINE
1848             size_type index1 () const {
1849                 return it_.index2 ();
1850             }
1851             BOOST_UBLAS_INLINE
1852             size_type index2 () const {
1853                 return it_.index1 ();
1854             }
1855
1856             // Assignment
1857             BOOST_UBLAS_INLINE
1858             const_iterator2 &operator = (const const_iterator2 &it) {
1859                 container_const_reference<self_type>::assign (&it ());
1860                 it_ = it.it_;
1861                 return *this;
1862             }
1863
1864             // Comparison
1865             BOOST_UBLAS_INLINE
1866             bool operator == (const const_iterator2 &it) const {
1867                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1868                 return it_ == it.it_;
1869             }
1870             BOOST_UBLAS_INLINE
1871             bool operator < (const const_iterator2 &it) const {
1872                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
1873                 return it_ < it.it_;
1874             }
1875
1876         private:
1877             const_subiterator2_type it_;
1878         };
1879 #endif
1880
1881         BOOST_UBLAS_INLINE
1882         const_iterator2 begin2 () const {
1883             return find2 (0, 0, 0);
1884         }
1885         BOOST_UBLAS_INLINE
1886         const_iterator2 cbegin2 () const {
1887             return begin2 ();
1888         }
1889         BOOST_UBLAS_INLINE
1890         const_iterator2 end2 () const {
1891             return find2 (0, 0, size2 ());
1892         }
1893         BOOST_UBLAS_INLINE
1894         const_iterator2 cend2 () const {
1895             return end2 ();
1896         }
1897
1898         // Reverse iterators
1899
1900         BOOST_UBLAS_INLINE
1901         const_reverse_iterator1 rbegin1 () const {
1902             return const_reverse_iterator1 (end1 ());
1903         }
1904         BOOST_UBLAS_INLINE
1905         const_reverse_iterator1 crbegin1 () const {
1906             return rbegin1 ();
1907         }
1908         BOOST_UBLAS_INLINE
1909         const_reverse_iterator1 rend1 () const {
1910             return const_reverse_iterator1 (begin1 ());
1911         }
1912         BOOST_UBLAS_INLINE
1913         const_reverse_iterator1 crend1 () const {
1914             return rend1 ();
1915         }
1916
1917         BOOST_UBLAS_INLINE
1918         const_reverse_iterator2 rbegin2 () const {
1919             return const_reverse_iterator2 (end2 ());
1920         }
1921         BOOST_UBLAS_INLINE
1922         const_reverse_iterator2 crbegin2 () const {
1923             return rbegin2 ();
1924         }
1925         BOOST_UBLAS_INLINE
1926         const_reverse_iterator2 rend2 () const {
1927             return const_reverse_iterator2 (begin2 ());
1928         }
1929         BOOST_UBLAS_INLINE
1930         const_reverse_iterator2 crend2 () const {
1931             return rend2 ();
1932         }
1933
1934     private:
1935         expression_closure_type e_;
1936     };
1937
1938     template<class E, class F>
1939     struct matrix_unary2_traits {
1940         typedef matrix_unary2<E, F> expression_type;
1941 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
1942         typedef expression_type result_type;
1943 #else
1944         typedef typename E::matrix_temporary_type result_type;
1945 #endif
1946     };
1947
1948     // (trans m) [i] [j] = m [j] [i]
1949     template<class E>
1950     BOOST_UBLAS_INLINE
1951     typename matrix_unary2_traits<const E, scalar_identity<typename E::value_type> >::result_type
1952     trans (const matrix_expression<E> &e) {
1953         typedef typename matrix_unary2_traits<const E, scalar_identity<typename E::value_type> >::expression_type expression_type;
1954         return expression_type (e ());
1955     }
1956     template<class E>
1957     BOOST_UBLAS_INLINE
1958     typename matrix_unary2_traits<E, scalar_identity<typename E::value_type> >::result_type
1959     trans (matrix_expression<E> &e) {
1960         typedef typename matrix_unary2_traits<E, scalar_identity<typename E::value_type> >::expression_type expression_type;
1961         return expression_type (e ());
1962     }
1963
1964     // (herm m) [i] [j] = conj (m [j] [i])
1965     template<class E>
1966     BOOST_UBLAS_INLINE
1967     typename matrix_unary2_traits<E, scalar_conj<typename E::value_type> >::result_type
1968     herm (const matrix_expression<E> &e) {
1969         typedef typename matrix_unary2_traits<E, scalar_conj<typename E::value_type> >::expression_type expression_type;
1970         return expression_type (e ());
1971     }
1972
1973     template<class E1, class E2, class F>
1974     class matrix_binary:
1975         public matrix_expression<matrix_binary<E1, E2, F> > {
1976
1977         typedef E1 expression1_type;
1978         typedef E2 expression2_type;
1979         typedef F functor_type;
1980     public:
1981         typedef typename E1::const_closure_type expression1_closure_type;
1982         typedef typename E2::const_closure_type expression2_closure_type;
1983     private:
1984         typedef matrix_binary<E1, E2, F> self_type;
1985     public:
1986 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
1987         using matrix_expression<self_type>::operator ();
1988 #endif
1989         typedef typename promote_traits<typename E1::size_type, typename E2::size_type>::promote_type size_type;
1990         typedef typename promote_traits<typename E1::difference_type, typename E2::difference_type>::promote_type difference_type;
1991         typedef typename F::result_type value_type;
1992         typedef value_type const_reference;
1993         typedef const_reference reference;
1994         typedef const self_type const_closure_type;
1995         typedef const_closure_type closure_type;
1996         typedef unknown_orientation_tag orientation_category;
1997         typedef unknown_storage_tag storage_category;
1998
1999         // Construction and destruction
2000         BOOST_UBLAS_INLINE
2001         matrix_binary (const E1 &e1, const E2 &e2):
2002             e1_ (e1), e2_ (e2) {}
2003
2004         // Accessors
2005         BOOST_UBLAS_INLINE
2006         size_type size1 () const {
2007             return BOOST_UBLAS_SAME (e1_.size1 (), e2_.size1 ());
2008         }
2009         BOOST_UBLAS_INLINE
2010         size_type size2 () const {
2011             return BOOST_UBLAS_SAME (e1_.size2 (), e2_.size2 ());
2012         }
2013
2014     public:
2015         // Expression accessors
2016         BOOST_UBLAS_INLINE
2017         const expression1_closure_type &expression1 () const {
2018             return e1_;
2019         }
2020         BOOST_UBLAS_INLINE
2021         const expression2_closure_type &expression2 () const {
2022             return e2_;
2023         }
2024
2025     public:
2026         // Element access
2027         BOOST_UBLAS_INLINE
2028         const_reference operator () (size_type i, size_type j) const {
2029             return functor_type::apply (e1_ (i, j), e2_ (i, j));
2030         }
2031
2032         // Closure comparison
2033         BOOST_UBLAS_INLINE
2034         bool same_closure (const matrix_binary &mb) const {
2035             return (*this).expression1 ().same_closure (mb.expression1 ()) &&
2036                    (*this).expression2 ().same_closure (mb.expression2 ());
2037         }
2038
2039         // Iterator types
2040     private:
2041         typedef typename E1::const_iterator1 const_iterator11_type;
2042         typedef typename E1::const_iterator2 const_iterator12_type;
2043         typedef typename E2::const_iterator1 const_iterator21_type;
2044         typedef typename E2::const_iterator2 const_iterator22_type;
2045         typedef const value_type *const_pointer;
2046
2047     public:
2048 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
2049         typedef typename iterator_restrict_traits<typename const_iterator11_type::iterator_category,
2050                                                   typename const_iterator21_type::iterator_category>::iterator_category iterator_category1;
2051         typedef indexed_const_iterator1<const_closure_type, iterator_category1> const_iterator1;
2052         typedef const_iterator1 iterator1;
2053         typedef typename iterator_restrict_traits<typename const_iterator12_type::iterator_category,
2054                                                   typename const_iterator22_type::iterator_category>::iterator_category iterator_category2;
2055         typedef indexed_const_iterator2<const_closure_type, iterator_category2> const_iterator2;
2056         typedef const_iterator2 iterator2;
2057 #else
2058         class const_iterator1;
2059         typedef const_iterator1 iterator1;
2060         class const_iterator2;
2061         typedef const_iterator2 iterator2;
2062 #endif
2063         typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
2064         typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
2065
2066         // Element lookup
2067         BOOST_UBLAS_INLINE
2068         const_iterator1 find1 (int rank, size_type i, size_type j) const {
2069             const_iterator11_type it11 (e1_.find1 (rank, i, j));
2070             const_iterator11_type it11_end (e1_.find1 (rank, size1 (), j));
2071             const_iterator21_type it21 (e2_.find1 (rank, i, j));
2072             const_iterator21_type it21_end (e2_.find1 (rank, size1 (), j));
2073             BOOST_UBLAS_CHECK (rank == 0 || it11 == it11_end || it11.index2 () == j, internal_logic ())
2074             BOOST_UBLAS_CHECK (rank == 0 || it21 == it21_end || it21.index2 () == j, internal_logic ())
2075             i = (std::min) (it11 != it11_end ? it11.index1 () : size1 (),
2076                           it21 != it21_end ? it21.index1 () : size1 ());
2077 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
2078             return const_iterator1 (*this, i, j);
2079 #else
2080             return const_iterator1 (*this, i, j, it11, it11_end, it21, it21_end);
2081 #endif
2082         }
2083         BOOST_UBLAS_INLINE
2084         const_iterator2 find2 (int rank, size_type i, size_type j) const {
2085             const_iterator12_type it12 (e1_.find2 (rank, i, j));
2086             const_iterator12_type it12_end (e1_.find2 (rank, i, size2 ()));
2087             const_iterator22_type it22 (e2_.find2 (rank, i, j));
2088             const_iterator22_type it22_end (e2_.find2 (rank, i, size2 ()));
2089             BOOST_UBLAS_CHECK (rank == 0 || it12 == it12_end || it12.index1 () == i, internal_logic ())
2090             BOOST_UBLAS_CHECK (rank == 0 || it22 == it22_end || it22.index1 () == i, internal_logic ())
2091             j = (std::min) (it12 != it12_end ? it12.index2 () : size2 (),
2092                           it22 != it22_end ? it22.index2 () : size2 ());
2093 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
2094             return const_iterator2 (*this, i, j);
2095 #else
2096             return const_iterator2 (*this, i, j, it12, it12_end, it22, it22_end);
2097 #endif
2098         }
2099
2100         // Iterators enhance the iterators of the referenced expression
2101         // with the binary functor.
2102
2103 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
2104         class const_iterator1:
2105             public container_const_reference<matrix_binary>,
2106             public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
2107                                                                           typename E2::const_iterator1::iterator_category>::iterator_category>::template
2108                 iterator_base<const_iterator1, value_type>::type {
2109         public:
2110             typedef typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
2111                                                       typename E2::const_iterator1::iterator_category>::iterator_category iterator_category;
2112             typedef typename matrix_binary::difference_type difference_type;
2113             typedef typename matrix_binary::value_type value_type;
2114             typedef typename matrix_binary::const_reference reference;
2115             typedef typename matrix_binary::const_pointer pointer;
2116
2117             typedef const_iterator2 dual_iterator_type;
2118             typedef const_reverse_iterator2 dual_reverse_iterator_type;
2119
2120             // Construction and destruction
2121             BOOST_UBLAS_INLINE
2122             const_iterator1 ():
2123                 container_const_reference<self_type> (), i_ (), j_ (), it1_ (), it1_end_ (), it2_ (), it2_end_ () {}
2124             BOOST_UBLAS_INLINE
2125             const_iterator1 (const self_type &mb, size_type i, size_type j,
2126                              const const_iterator11_type &it1, const const_iterator11_type &it1_end,
2127                              const const_iterator21_type &it2, const const_iterator21_type &it2_end):
2128                 container_const_reference<self_type> (mb), i_ (i), j_ (j), it1_ (it1), it1_end_ (it1_end), it2_ (it2), it2_end_ (it2_end) {}
2129
2130         private:
2131             // Dense specializations
2132             BOOST_UBLAS_INLINE
2133             void increment (dense_random_access_iterator_tag) {
2134                 ++ i_; ++ it1_; ++ it2_;
2135             }
2136             BOOST_UBLAS_INLINE
2137             void decrement (dense_random_access_iterator_tag) {
2138                 -- i_; -- it1_; -- it2_;
2139             }
2140             BOOST_UBLAS_INLINE
2141             void increment (dense_random_access_iterator_tag, difference_type n) {
2142                  i_ += n; it1_ += n; it2_ += n;
2143             }
2144             BOOST_UBLAS_INLINE
2145             void decrement (dense_random_access_iterator_tag, difference_type n) {
2146                 i_ -= n; it1_ -= n; it2_ -= n;
2147             }
2148             BOOST_UBLAS_INLINE
2149             value_type dereference (dense_random_access_iterator_tag) const {
2150                 return functor_type::apply (*it1_, *it2_);
2151             }
2152
2153             // Packed specializations
2154             BOOST_UBLAS_INLINE
2155             void increment (packed_random_access_iterator_tag) {
2156                 if (it1_ != it1_end_)
2157                     if (it1_.index1 () <= i_)
2158                         ++ it1_;
2159                 if (it2_ != it2_end_)
2160                     if (it2_.index1 () <= i_)
2161                         ++ it2_;
2162                 ++ i_;
2163             }
2164             BOOST_UBLAS_INLINE
2165             void decrement (packed_random_access_iterator_tag) {
2166                 if (it1_ != it1_end_)
2167                     if (i_ <= it1_.index1 ())
2168                         -- it1_;
2169                 if (it2_ != it2_end_)
2170                     if (i_ <= it2_.index1 ())
2171                         -- it2_;
2172                 -- i_;
2173             }
2174             BOOST_UBLAS_INLINE
2175             void increment (packed_random_access_iterator_tag, difference_type n) {
2176                 while (n > 0) {
2177                     increment (packed_random_access_iterator_tag ());
2178                     --n;
2179                 }
2180                 while (n < 0) {
2181                     decrement (packed_random_access_iterator_tag ());
2182                     ++n;
2183                 }
2184             }
2185             BOOST_UBLAS_INLINE
2186             void decrement (packed_random_access_iterator_tag, difference_type n) {
2187                 while (n > 0) {
2188                     decrement (packed_random_access_iterator_tag ());
2189                     --n;
2190                 }
2191                 while (n < 0) {
2192                     increment (packed_random_access_iterator_tag ());
2193                     ++n;
2194                 }
2195             }
2196             BOOST_UBLAS_INLINE
2197             value_type dereference (packed_random_access_iterator_tag) const {
2198                 value_type t1 = value_type/*zero*/();
2199                 if (it1_ != it1_end_) {
2200                     BOOST_UBLAS_CHECK (it1_.index2 () == j_, internal_logic ());
2201                     if (it1_.index1 () == i_)
2202                         t1 = *it1_;
2203                 }
2204                 value_type t2 = value_type/*zero*/();
2205                 if (it2_ != it2_end_) {
2206                     BOOST_UBLAS_CHECK (it2_.index2 () == j_, internal_logic ());
2207                     if (it2_.index1 () == i_)
2208                         t2 = *it2_;
2209                 }
2210                 return functor_type::apply (t1, t2);
2211             }
2212
2213             // Sparse specializations
2214             BOOST_UBLAS_INLINE
2215             void increment (sparse_bidirectional_iterator_tag) {
2216                 size_type index1 = (*this) ().size1 ();
2217                 if (it1_ != it1_end_) {
2218                     if (it1_.index1 () <= i_)
2219                         ++ it1_;
2220                     if (it1_ != it1_end_)
2221                         index1 = it1_.index1 ();
2222                 }
2223                 size_type index2 = (*this) ().size1 ();
2224                 if (it2_ != it2_end_)
2225                     if (it2_.index1 () <= i_)
2226                         ++ it2_;
2227                     if (it2_ != it2_end_) {
2228                         index2 = it2_.index1 ();
2229                 }
2230                 i_ = (std::min) (index1, index2);
2231             }
2232             BOOST_UBLAS_INLINE
2233             void decrement (sparse_bidirectional_iterator_tag) {
2234                 size_type index1 = (*this) ().size1 ();
2235                 if (it1_ != it1_end_) {
2236                     if (i_ <= it1_.index1 ())
2237                         -- it1_;
2238                     if (it1_ != it1_end_)
2239                         index1 = it1_.index1 ();
2240                 }
2241                 size_type index2 = (*this) ().size1 ();
2242                 if (it2_ != it2_end_) {
2243                     if (i_ <= it2_.index1 ())
2244                         -- it2_;
2245                     if (it2_ != it2_end_)
2246                         index2 = it2_.index1 ();
2247                 }
2248                 i_ = (std::max) (index1, index2);
2249             }
2250             BOOST_UBLAS_INLINE
2251             void increment (sparse_bidirectional_iterator_tag, difference_type n) {
2252                 while (n > 0) {
2253                     increment (sparse_bidirectional_iterator_tag ());
2254                     --n;
2255                 }
2256                 while (n < 0) {
2257                     decrement (sparse_bidirectional_iterator_tag ());
2258                     ++n;
2259                 }
2260             }
2261             BOOST_UBLAS_INLINE
2262             void decrement (sparse_bidirectional_iterator_tag, difference_type n) {
2263                 while (n > 0) {
2264                     decrement (sparse_bidirectional_iterator_tag ());
2265                     --n;
2266                 }
2267                 while (n < 0) {
2268                     increment (sparse_bidirectional_iterator_tag ());
2269                     ++n;
2270                 }
2271             }
2272             BOOST_UBLAS_INLINE
2273             value_type dereference (sparse_bidirectional_iterator_tag) const {
2274                 value_type t1 = value_type/*zero*/();
2275                 if (it1_ != it1_end_) {
2276                     BOOST_UBLAS_CHECK (it1_.index2 () == j_, internal_logic ());
2277                     if (it1_.index1 () == i_)
2278                         t1 = *it1_;
2279                 }
2280                 value_type t2 = value_type/*zero*/();
2281                 if (it2_ != it2_end_) {
2282                     BOOST_UBLAS_CHECK (it2_.index2 () == j_, internal_logic ());
2283                     if (it2_.index1 () == i_)
2284                         t2 = *it2_;
2285                 }
2286                 return functor_type::apply (t1, t2);
2287             }
2288
2289         public:
2290             // Arithmetic
2291             BOOST_UBLAS_INLINE
2292             const_iterator1 &operator ++ () {
2293                 increment (iterator_category ());
2294                 return *this;
2295             }
2296             BOOST_UBLAS_INLINE
2297             const_iterator1 &operator -- () {
2298                 decrement (iterator_category ());
2299                 return *this;
2300             }
2301             BOOST_UBLAS_INLINE
2302             const_iterator1 &operator += (difference_type n) {
2303                 increment (iterator_category (), n);
2304                 return *this;
2305             }
2306             BOOST_UBLAS_INLINE
2307             const_iterator1 &operator -= (difference_type n) {
2308                 decrement (iterator_category (), n);
2309                 return *this;
2310             }
2311             BOOST_UBLAS_INLINE
2312             difference_type operator - (const const_iterator1 &it) const {
2313                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2314                 BOOST_UBLAS_CHECK (index2 () == it.index2 (), external_logic ());
2315                 return index1 () - it.index1 ();
2316             }
2317
2318             // Dereference
2319             BOOST_UBLAS_INLINE
2320             const_reference operator * () const {
2321                 return dereference (iterator_category ());
2322             }
2323             BOOST_UBLAS_INLINE
2324             const_reference operator [] (difference_type n) const {
2325                 return *(*this + n);
2326             }
2327
2328 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
2329             BOOST_UBLAS_INLINE
2330 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2331             typename self_type::
2332 #endif
2333             const_iterator2 begin () const {
2334                 return (*this) ().find2 (1, index1 (), 0);
2335             }
2336             BOOST_UBLAS_INLINE
2337 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2338             typename self_type::
2339 #endif
2340             const_iterator2 cbegin () const {
2341                 return begin ();
2342             }
2343             BOOST_UBLAS_INLINE
2344 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2345             typename self_type::
2346 #endif
2347             const_iterator2 end () const {
2348                 return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
2349             }
2350             BOOST_UBLAS_INLINE
2351 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2352             typename self_type::
2353 #endif
2354             const_iterator2 cend () const {
2355                 return end ();
2356             }
2357             BOOST_UBLAS_INLINE
2358 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2359             typename self_type::
2360 #endif
2361             const_reverse_iterator2 rbegin () const {
2362                 return const_reverse_iterator2 (end ());
2363             }
2364             BOOST_UBLAS_INLINE
2365 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2366             typename self_type::
2367 #endif
2368             const_reverse_iterator2 crbegin () const {
2369                 return rbegin ();
2370             }
2371             BOOST_UBLAS_INLINE
2372 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2373             typename self_type::
2374 #endif
2375             const_reverse_iterator2 rend () const {
2376                 return const_reverse_iterator2 (begin ());
2377             }
2378             BOOST_UBLAS_INLINE
2379 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2380             typename self_type::
2381 #endif
2382             const_reverse_iterator2 crend () const {
2383                 return rend ();
2384             }
2385 #endif
2386
2387             // Indices
2388             BOOST_UBLAS_INLINE
2389             size_type index1 () const {
2390                 return i_;
2391             }
2392             BOOST_UBLAS_INLINE
2393             size_type index2 () const {
2394                 // if (it1_ != it1_end_ && it2_ != it2_end_)
2395                 //    return BOOST_UBLAS_SAME (it1_.index2 (), it2_.index2 ());
2396                 // else
2397                     return j_;
2398             }
2399
2400             // Assignment
2401             BOOST_UBLAS_INLINE
2402             const_iterator1 &operator = (const const_iterator1 &it) {
2403                 container_const_reference<self_type>::assign (&it ());
2404                 i_ = it.i_;
2405                 j_ = it.j_;
2406                 it1_ = it.it1_;
2407                 it1_end_ = it.it1_end_;
2408                 it2_ = it.it2_;
2409                 it2_end_ = it.it2_end_;
2410                 return *this;
2411             }
2412
2413             // Comparison
2414             BOOST_UBLAS_INLINE
2415             bool operator == (const const_iterator1 &it) const {
2416                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2417                 BOOST_UBLAS_CHECK (index2 () == it.index2 (), external_logic ());
2418                 return index1 () == it.index1 ();
2419             }
2420             BOOST_UBLAS_INLINE
2421             bool operator < (const const_iterator1 &it) const {
2422                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2423                 BOOST_UBLAS_CHECK (index2 () == it.index2 (), external_logic ());
2424                 return index1 () < it.index1 ();
2425             }
2426
2427         private:
2428             size_type i_;
2429             size_type j_;
2430             const_iterator11_type it1_;
2431             const_iterator11_type it1_end_;
2432             const_iterator21_type it2_;
2433             const_iterator21_type it2_end_;
2434         };
2435 #endif
2436
2437         BOOST_UBLAS_INLINE
2438         const_iterator1 begin1 () const {
2439             return find1 (0, 0, 0);
2440         }
2441         BOOST_UBLAS_INLINE
2442         const_iterator1 cbegin1 () const {
2443             return begin1 ();
2444         }
2445         BOOST_UBLAS_INLINE
2446         const_iterator1 end1 () const {
2447             return find1 (0, size1 (), 0);
2448         }
2449         BOOST_UBLAS_INLINE
2450         const_iterator1 cend1 () const {
2451             return end1 ();
2452         }
2453
2454 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
2455         class const_iterator2:
2456             public container_const_reference<matrix_binary>,
2457             public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator2::iterator_category,
2458                                                                           typename E2::const_iterator2::iterator_category>::iterator_category>::template
2459                 iterator_base<const_iterator2, value_type>::type {
2460         public:
2461             typedef typename iterator_restrict_traits<typename E1::const_iterator2::iterator_category,
2462                                                       typename E2::const_iterator2::iterator_category>::iterator_category iterator_category;
2463             typedef typename matrix_binary::difference_type difference_type;
2464             typedef typename matrix_binary::value_type value_type;
2465             typedef typename matrix_binary::const_reference reference;
2466             typedef typename matrix_binary::const_pointer pointer;
2467
2468             typedef const_iterator1 dual_iterator_type;
2469             typedef const_reverse_iterator1 dual_reverse_iterator_type;
2470
2471             // Construction and destruction
2472             BOOST_UBLAS_INLINE
2473             const_iterator2 ():
2474                 container_const_reference<self_type> (), i_ (), j_ (), it1_ (), it1_end_ (), it2_ (), it2_end_ () {}
2475             BOOST_UBLAS_INLINE
2476             const_iterator2 (const self_type &mb, size_type i, size_type j,
2477                              const const_iterator12_type &it1, const const_iterator12_type &it1_end,
2478                              const const_iterator22_type &it2, const const_iterator22_type &it2_end):
2479                 container_const_reference<self_type> (mb), i_ (i), j_ (j), it1_ (it1), it1_end_ (it1_end), it2_ (it2), it2_end_ (it2_end) {}
2480
2481         private:
2482             // Dense access specializations
2483             BOOST_UBLAS_INLINE
2484             void increment (dense_random_access_iterator_tag) {
2485                 ++ j_; ++ it1_; ++ it2_;
2486             }
2487             BOOST_UBLAS_INLINE
2488             void decrement (dense_random_access_iterator_tag) {
2489                 -- j_; -- it1_; -- it2_;
2490             }
2491             BOOST_UBLAS_INLINE
2492             void increment (dense_random_access_iterator_tag, difference_type n) {
2493                 j_ += n; it1_ += n; it2_ += n;
2494             }
2495             BOOST_UBLAS_INLINE
2496             void decrement (dense_random_access_iterator_tag, difference_type n) {
2497                 j_ -= n; it1_ -= n; it2_ -= n;
2498             }
2499             BOOST_UBLAS_INLINE
2500             value_type dereference (dense_random_access_iterator_tag) const {
2501                 return functor_type::apply (*it1_, *it2_);
2502             }
2503
2504             // Packed specializations
2505             BOOST_UBLAS_INLINE
2506             void increment (packed_random_access_iterator_tag) {
2507                 if (it1_ != it1_end_)
2508                     if (it1_.index2 () <= j_)
2509                         ++ it1_;
2510                 if (it2_ != it2_end_)
2511                     if (it2_.index2 () <= j_)
2512                         ++ it2_;
2513                 ++ j_;
2514             }
2515             BOOST_UBLAS_INLINE
2516             void decrement (packed_random_access_iterator_tag) {
2517                 if (it1_ != it1_end_)
2518                     if (j_ <= it1_.index2 ())
2519                         -- it1_;
2520                 if (it2_ != it2_end_)
2521                     if (j_ <= it2_.index2 ())
2522                         -- it2_;
2523                 -- j_;
2524             }
2525             BOOST_UBLAS_INLINE
2526             void increment (packed_random_access_iterator_tag, difference_type n) {
2527                 while (n > 0) {
2528                     increment (packed_random_access_iterator_tag ());
2529                     --n;
2530                 }
2531                 while (n < 0) {
2532                     decrement (packed_random_access_iterator_tag ());
2533                     ++n;
2534                 }
2535             }
2536             BOOST_UBLAS_INLINE
2537             void decrement (packed_random_access_iterator_tag, difference_type n) {
2538                 while (n > 0) {
2539                     decrement (packed_random_access_iterator_tag ());
2540                     --n;
2541                 }
2542                 while (n < 0) {
2543                     increment (packed_random_access_iterator_tag ());
2544                     ++n;
2545                 }
2546             }
2547             BOOST_UBLAS_INLINE
2548             value_type dereference (packed_random_access_iterator_tag) const {
2549                 value_type t1 = value_type/*zero*/();
2550                 if (it1_ != it1_end_) {
2551                     BOOST_UBLAS_CHECK (it1_.index1 () == i_, internal_logic ());
2552                     if (it1_.index2 () == j_)
2553                         t1 = *it1_;
2554                 }
2555                 value_type t2 = value_type/*zero*/();
2556                 if (it2_ != it2_end_) {
2557                     BOOST_UBLAS_CHECK (it2_.index1 () == i_, internal_logic ());
2558                     if (it2_.index2 () == j_)
2559                         t2 = *it2_;
2560                 }
2561                 return functor_type::apply (t1, t2);
2562             }
2563
2564             // Sparse specializations
2565             BOOST_UBLAS_INLINE
2566             void increment (sparse_bidirectional_iterator_tag) {
2567                 size_type index1 = (*this) ().size2 ();
2568                 if (it1_ != it1_end_) {
2569                     if (it1_.index2 () <= j_)
2570                         ++ it1_;
2571                     if (it1_ != it1_end_)
2572                         index1 = it1_.index2 ();
2573                 }
2574                 size_type index2 = (*this) ().size2 ();
2575                 if (it2_ != it2_end_) {
2576                     if (it2_.index2 () <= j_)
2577                         ++ it2_;
2578                     if (it2_ != it2_end_)
2579                         index2 = it2_.index2 ();
2580                 }
2581                 j_ = (std::min) (index1, index2);
2582             }
2583             BOOST_UBLAS_INLINE
2584             void decrement (sparse_bidirectional_iterator_tag) {
2585                 size_type index1 = (*this) ().size2 ();
2586                 if (it1_ != it1_end_) {
2587                     if (j_ <= it1_.index2 ())
2588                         -- it1_;
2589                     if (it1_ != it1_end_)
2590                         index1 = it1_.index2 ();
2591                 }
2592                 size_type index2 = (*this) ().size2 ();
2593                 if (it2_ != it2_end_) {
2594                     if (j_ <= it2_.index2 ())
2595                         -- it2_;
2596                     if (it2_ != it2_end_)
2597                         index2 = it2_.index2 ();
2598                 }
2599                 j_ = (std::max) (index1, index2);
2600             }
2601             BOOST_UBLAS_INLINE
2602             void increment (sparse_bidirectional_iterator_tag, difference_type n) {
2603                 while (n > 0) {
2604                     increment (sparse_bidirectional_iterator_tag ());
2605                     --n;
2606                 }
2607                 while (n < 0) {
2608                     decrement (sparse_bidirectional_iterator_tag ());
2609                     ++n;
2610                 }
2611             }
2612             BOOST_UBLAS_INLINE
2613             void decrement (sparse_bidirectional_iterator_tag, difference_type n) {
2614                 while (n > 0) {
2615                     decrement (sparse_bidirectional_iterator_tag ());
2616                     --n;
2617                 }
2618                 while (n < 0) {
2619                     increment (sparse_bidirectional_iterator_tag ());
2620                     ++n;
2621                 }
2622             }
2623             BOOST_UBLAS_INLINE
2624             value_type dereference (sparse_bidirectional_iterator_tag) const {
2625                 value_type t1 = value_type/*zero*/();
2626                 if (it1_ != it1_end_) {
2627                     BOOST_UBLAS_CHECK (it1_.index1 () == i_, internal_logic ());
2628                     if (it1_.index2 () == j_)
2629                         t1 = *it1_;
2630                 }
2631                 value_type t2 = value_type/*zero*/();
2632                 if (it2_ != it2_end_) {
2633                     BOOST_UBLAS_CHECK (it2_.index1 () == i_, internal_logic ());
2634                     if (it2_.index2 () == j_)
2635                         t2 = *it2_;
2636                 }
2637                 return functor_type::apply (t1, t2);
2638             }
2639
2640         public:
2641             // Arithmetic
2642             BOOST_UBLAS_INLINE
2643             const_iterator2 &operator ++ () {
2644                 increment (iterator_category ());
2645                 return *this;
2646             }
2647             BOOST_UBLAS_INLINE
2648             const_iterator2 &operator -- () {
2649                 decrement (iterator_category ());
2650                 return *this;
2651             }
2652             BOOST_UBLAS_INLINE
2653             const_iterator2 &operator += (difference_type n) {
2654                 increment (iterator_category (), n);
2655                 return *this;
2656             }
2657             BOOST_UBLAS_INLINE
2658             const_iterator2 &operator -= (difference_type n) {
2659                 decrement (iterator_category (), n);
2660                 return *this;
2661             }
2662             BOOST_UBLAS_INLINE
2663             difference_type operator - (const const_iterator2 &it) const {
2664                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2665                 BOOST_UBLAS_CHECK (index1 () == it.index1 (), external_logic ());
2666                 return index2 () - it.index2 ();
2667             }
2668
2669             // Dereference
2670             BOOST_UBLAS_INLINE
2671             const_reference operator * () const {
2672                 return dereference (iterator_category ());
2673             }
2674             BOOST_UBLAS_INLINE
2675             const_reference operator [] (difference_type n) const {
2676                 return *(*this + n);
2677             }
2678
2679 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
2680             BOOST_UBLAS_INLINE
2681 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2682             typename self_type::
2683 #endif
2684             const_iterator1 begin () const {
2685                 return (*this) ().find1 (1, 0, index2 ());
2686             }
2687             BOOST_UBLAS_INLINE
2688 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2689             typename self_type::
2690 #endif
2691             const_iterator1 cbegin () const {
2692                 return begin ();
2693             }
2694             BOOST_UBLAS_INLINE
2695 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2696             typename self_type::
2697 #endif
2698             const_iterator1 end () const {
2699                 return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
2700             }
2701             BOOST_UBLAS_INLINE
2702 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2703             typename self_type::
2704 #endif
2705             const_iterator1 cend () const {
2706                 return end ();
2707             }
2708             BOOST_UBLAS_INLINE
2709 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2710             typename self_type::
2711 #endif
2712             const_reverse_iterator1 rbegin () const {
2713                 return const_reverse_iterator1 (end ());
2714             }
2715             BOOST_UBLAS_INLINE
2716 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2717             typename self_type::
2718 #endif
2719             const_reverse_iterator1 crbegin () const {
2720                 return rbegin ();
2721             }
2722             BOOST_UBLAS_INLINE
2723 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2724             typename self_type::
2725 #endif
2726             const_reverse_iterator1 rend () const {
2727                 return const_reverse_iterator1 (begin ());
2728             }
2729             BOOST_UBLAS_INLINE
2730 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2731             typename self_type::
2732 #endif
2733             const_reverse_iterator1 crend () const {
2734                 return rend ();
2735             }
2736 #endif
2737
2738             // Indices
2739             BOOST_UBLAS_INLINE
2740             size_type index1 () const {
2741                 // if (it1_ != it1_end_ && it2_ != it2_end_)
2742                 //    return BOOST_UBLAS_SAME (it1_.index1 (), it2_.index1 ());
2743                 // else
2744                     return i_;
2745             }
2746             BOOST_UBLAS_INLINE
2747             size_type index2 () const {
2748                 return j_;
2749             }
2750
2751             // Assignment
2752             BOOST_UBLAS_INLINE
2753             const_iterator2 &operator = (const const_iterator2 &it) {
2754                 container_const_reference<self_type>::assign (&it ());
2755                 i_ = it.i_;
2756                 j_ = it.j_;
2757                 it1_ = it.it1_;
2758                 it1_end_ = it.it1_end_;
2759                 it2_ = it.it2_;
2760                 it2_end_ = it.it2_end_;
2761                 return *this;
2762             }
2763
2764             // Comparison
2765             BOOST_UBLAS_INLINE
2766             bool operator == (const const_iterator2 &it) const {
2767                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2768                 BOOST_UBLAS_CHECK (index1 () == it.index1 (), external_logic ());
2769                 return index2 () == it.index2 ();
2770             }
2771             BOOST_UBLAS_INLINE
2772             bool operator < (const const_iterator2 &it) const {
2773                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
2774                 BOOST_UBLAS_CHECK (index1 () == it.index1 (), external_logic ());
2775                 return index2 () < it.index2 ();
2776             }
2777
2778         private:
2779             size_type i_;
2780             size_type j_;
2781             const_iterator12_type it1_;
2782             const_iterator12_type it1_end_;
2783             const_iterator22_type it2_;
2784             const_iterator22_type it2_end_;
2785         };
2786 #endif
2787
2788         BOOST_UBLAS_INLINE
2789         const_iterator2 begin2 () const {
2790             return find2 (0, 0, 0);
2791         }
2792         BOOST_UBLAS_INLINE
2793         const_iterator2 cbegin2 () const {
2794             return begin2 ();
2795         }
2796         BOOST_UBLAS_INLINE
2797         const_iterator2 end2 () const {
2798             return find2 (0, 0, size2 ());
2799         }
2800         BOOST_UBLAS_INLINE
2801         const_iterator2 cend2 () const {
2802             return end2 ();
2803         }
2804
2805         // Reverse iterators
2806
2807         BOOST_UBLAS_INLINE
2808         const_reverse_iterator1 rbegin1 () const {
2809             return const_reverse_iterator1 (end1 ());
2810         }
2811         BOOST_UBLAS_INLINE
2812         const_reverse_iterator1 crbegin1 () const {
2813             return rbegin1 ();
2814         }
2815         BOOST_UBLAS_INLINE
2816         const_reverse_iterator1 rend1 () const {
2817             return const_reverse_iterator1 (begin1 ());
2818         }
2819         BOOST_UBLAS_INLINE
2820         const_reverse_iterator1 crend1 () const {
2821             return rend1 ();
2822         }
2823
2824         BOOST_UBLAS_INLINE
2825         const_reverse_iterator2 rbegin2 () const {
2826             return const_reverse_iterator2 (end2 ());
2827         }
2828         BOOST_UBLAS_INLINE
2829         const_reverse_iterator2 crbegin2 () const {
2830             return rbegin2 ();
2831         }
2832         BOOST_UBLAS_INLINE
2833         const_reverse_iterator2 rend2 () const {
2834             return const_reverse_iterator2 (begin2 ());
2835         }
2836         BOOST_UBLAS_INLINE
2837         const_reverse_iterator2 crend2 () const {
2838             return rend2 ();
2839         }
2840
2841     private:
2842         expression1_closure_type e1_;
2843         expression2_closure_type e2_;
2844     };
2845
2846     template<class E1, class E2, class F>
2847     struct matrix_binary_traits {
2848         typedef matrix_binary<E1, E2, F> expression_type;
2849 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
2850         typedef expression_type result_type;
2851 #else
2852         typedef typename E1::matrix_temporary_type result_type;
2853 #endif
2854     };
2855
2856     // (m1 + m2) [i] [j] = m1 [i] [j] + m2 [i] [j]
2857     template<class E1, class E2>
2858     BOOST_UBLAS_INLINE
2859     typename matrix_binary_traits<E1, E2, scalar_plus<typename E1::value_type,
2860                                                       typename E2::value_type> >::result_type
2861     operator + (const matrix_expression<E1> &e1,
2862                 const matrix_expression<E2> &e2) {
2863         typedef typename matrix_binary_traits<E1, E2, scalar_plus<typename E1::value_type,
2864                                                                   typename E2::value_type> >::expression_type expression_type;
2865         return expression_type (e1 (), e2 ());
2866     }
2867
2868     // (m1 - m2) [i] [j] = m1 [i] [j] - m2 [i] [j]
2869     template<class E1, class E2>
2870     BOOST_UBLAS_INLINE
2871     typename matrix_binary_traits<E1, E2, scalar_minus<typename E1::value_type,
2872                                                        typename E2::value_type> >::result_type
2873     operator - (const matrix_expression<E1> &e1,
2874                 const matrix_expression<E2> &e2) {
2875         typedef typename matrix_binary_traits<E1, E2, scalar_minus<typename E1::value_type,
2876                                                                    typename E2::value_type> >::expression_type expression_type;
2877         return expression_type (e1 (), e2 ());
2878     }
2879
2880     // (m1 * m2) [i] [j] = m1 [i] [j] * m2 [i] [j]
2881     template<class E1, class E2>
2882     BOOST_UBLAS_INLINE
2883     typename matrix_binary_traits<E1, E2, scalar_multiplies<typename E1::value_type,
2884                                                             typename E2::value_type> >::result_type
2885     element_prod (const matrix_expression<E1> &e1,
2886                   const matrix_expression<E2> &e2) {
2887         typedef typename matrix_binary_traits<E1, E2, scalar_multiplies<typename E1::value_type,
2888                                                                         typename E2::value_type> >::expression_type expression_type;
2889         return expression_type (e1 (), e2 ());
2890     }
2891
2892     // (m1 / m2) [i] [j] = m1 [i] [j] / m2 [i] [j]
2893     template<class E1, class E2>
2894     BOOST_UBLAS_INLINE
2895     typename matrix_binary_traits<E1, E2, scalar_divides<typename E1::value_type,
2896                                                          typename E2::value_type> >::result_type
2897     element_div (const matrix_expression<E1> &e1,
2898                  const matrix_expression<E2> &e2) {
2899         typedef typename matrix_binary_traits<E1, E2, scalar_divides<typename E1::value_type,
2900                                                                      typename E2::value_type> >::expression_type expression_type;
2901         return expression_type (e1 (), e2 ());
2902     }
2903
2904     template<class E1, class E2, class F>
2905     class matrix_binary_scalar1:
2906         public matrix_expression<matrix_binary_scalar1<E1, E2, F> > {
2907
2908         typedef E1 expression1_type;
2909         typedef E2 expression2_type;
2910         typedef F functor_type;
2911         typedef const E1& expression1_closure_type;
2912         typedef typename E2::const_closure_type expression2_closure_type;
2913         typedef matrix_binary_scalar1<E1, E2, F> self_type;
2914     public:
2915 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
2916         using matrix_expression<self_type>::operator ();
2917 #endif
2918         typedef typename E2::size_type size_type;
2919         typedef typename E2::difference_type difference_type;
2920         typedef typename F::result_type value_type;
2921         typedef value_type const_reference;
2922         typedef const_reference reference;
2923         typedef const self_type const_closure_type;
2924         typedef const_closure_type closure_type;
2925         typedef typename E2::orientation_category orientation_category;
2926         typedef unknown_storage_tag storage_category;
2927
2928         // Construction and destruction
2929         BOOST_UBLAS_INLINE
2930         matrix_binary_scalar1 (const expression1_type &e1, const expression2_type &e2):
2931             e1_ (e1), e2_ (e2) {}
2932
2933         // Accessors
2934         BOOST_UBLAS_INLINE
2935         size_type size1 () const {
2936             return e2_.size1 ();
2937         }
2938         BOOST_UBLAS_INLINE
2939         size_type size2 () const {
2940             return e2_.size2 ();
2941         }
2942
2943     public:
2944         // Element access
2945         BOOST_UBLAS_INLINE
2946         const_reference operator () (size_type i, size_type j) const {
2947             return functor_type::apply (expression1_type (e1_), e2_ (i, j));
2948         }
2949
2950         // Closure comparison
2951         BOOST_UBLAS_INLINE
2952         bool same_closure (const matrix_binary_scalar1 &mbs1) const {
2953             return &e1_ == &(mbs1.e1_) &&
2954                    (*this).e2_.same_closure (mbs1.e2_);
2955         }
2956
2957         // Iterator types
2958     private:
2959         typedef expression1_type const_subiterator1_type;
2960         typedef typename E2::const_iterator1 const_iterator21_type;
2961         typedef typename E2::const_iterator2 const_iterator22_type;
2962         typedef const value_type *const_pointer;
2963
2964     public:
2965 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
2966         typedef indexed_const_iterator1<const_closure_type, typename const_iterator21_type::iterator_category> const_iterator1;
2967         typedef const_iterator1 iterator1;
2968         typedef indexed_const_iterator2<const_closure_type, typename const_iterator22_type::iterator_category> const_iterator2;
2969         typedef const_iterator2 iterator2;
2970 #else
2971         class const_iterator1;
2972         typedef const_iterator1 iterator1;
2973         class const_iterator2;
2974         typedef const_iterator2 iterator2;
2975 #endif
2976         typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
2977         typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
2978
2979         // Element lookup
2980         BOOST_UBLAS_INLINE
2981         const_iterator1 find1 (int rank, size_type i, size_type j) const {
2982             const_iterator21_type it21 (e2_.find1 (rank, i, j));
2983 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
2984             return const_iterator1 (*this, it21.index1 (), it21.index2 ());
2985 #else
2986             return const_iterator1 (*this, const_subiterator1_type (e1_), it21);
2987 #endif
2988         }
2989         BOOST_UBLAS_INLINE
2990         const_iterator2 find2 (int rank, size_type i, size_type j) const {
2991             const_iterator22_type it22 (e2_.find2 (rank, i, j));
2992 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
2993             return const_iterator2 (*this, it22.index1 (), it22.index2 ());
2994 #else
2995             return const_iterator2 (*this, const_subiterator1_type (e1_), it22);
2996 #endif
2997         }
2998
2999         // Iterators enhance the iterators of the referenced expression
3000         // with the binary functor.
3001
3002 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
3003         class const_iterator1:
3004             public container_const_reference<matrix_binary_scalar1>,
3005             public iterator_base_traits<typename E2::const_iterator1::iterator_category>::template
3006                 iterator_base<const_iterator1, value_type>::type {
3007         public:
3008             typedef typename E2::const_iterator1::iterator_category iterator_category;
3009             typedef typename matrix_binary_scalar1::difference_type difference_type;
3010             typedef typename matrix_binary_scalar1::value_type value_type;
3011             typedef typename matrix_binary_scalar1::const_reference reference;
3012             typedef typename matrix_binary_scalar1::const_pointer pointer;
3013
3014             typedef const_iterator2 dual_iterator_type;
3015             typedef const_reverse_iterator2 dual_reverse_iterator_type;
3016
3017             // Construction and destruction
3018             BOOST_UBLAS_INLINE
3019             const_iterator1 ():
3020                 container_const_reference<self_type> (), it1_ (), it2_ () {}
3021             BOOST_UBLAS_INLINE
3022             const_iterator1 (const self_type &mbs, const const_subiterator1_type &it1, const const_iterator21_type &it2):
3023                 container_const_reference<self_type> (mbs), it1_ (it1), it2_ (it2) {}
3024
3025             // Arithmetic
3026             BOOST_UBLAS_INLINE
3027             const_iterator1 &operator ++ () {
3028                 ++ it2_;
3029                 return *this;
3030             }
3031             BOOST_UBLAS_INLINE
3032             const_iterator1 &operator -- () {
3033                 -- it2_ ;
3034                 return *this;
3035             }
3036             BOOST_UBLAS_INLINE
3037             const_iterator1 &operator += (difference_type n) {
3038                 it2_ += n;
3039                 return *this;
3040             }
3041             BOOST_UBLAS_INLINE
3042             const_iterator1 &operator -= (difference_type n) {
3043                 it2_ -= n;
3044                 return *this;
3045             }
3046             BOOST_UBLAS_INLINE
3047             difference_type operator - (const const_iterator1 &it) const {
3048                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3049                 // FIXME we shouldn't compare floats
3050                 // BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
3051                 return it2_ - it.it2_;
3052             }
3053
3054             // Dereference
3055             BOOST_UBLAS_INLINE
3056             const_reference operator * () const {
3057                 return functor_type::apply (it1_, *it2_);
3058             }
3059             BOOST_UBLAS_INLINE
3060             const_reference operator [] (difference_type n) const {
3061                 return *(*this + n);
3062             }
3063
3064 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
3065             BOOST_UBLAS_INLINE
3066 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3067             typename self_type::
3068 #endif
3069             const_iterator2 begin () const {
3070                 return (*this) ().find2 (1, index1 (), 0);
3071             }
3072             BOOST_UBLAS_INLINE
3073 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3074             typename self_type::
3075 #endif
3076             const_iterator2 cbegin () const {
3077                 return begin ();
3078             }
3079             BOOST_UBLAS_INLINE
3080 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3081             typename self_type::
3082 #endif
3083             const_iterator2 end () const {
3084                 return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
3085             }
3086             BOOST_UBLAS_INLINE
3087 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3088             typename self_type::
3089 #endif
3090             const_iterator2 cend () const {
3091                 return end ();
3092             }
3093             BOOST_UBLAS_INLINE
3094 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3095             typename self_type::
3096 #endif
3097             const_reverse_iterator2 rbegin () const {
3098                 return const_reverse_iterator2 (end ());
3099             }
3100             BOOST_UBLAS_INLINE
3101 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3102             typename self_type::
3103 #endif
3104             const_reverse_iterator2 crbegin () const {
3105                 return rbegin ();
3106             }
3107             BOOST_UBLAS_INLINE
3108 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3109             typename self_type::
3110 #endif
3111             const_reverse_iterator2 rend () const {
3112                 return const_reverse_iterator2 (begin ());
3113             }
3114             BOOST_UBLAS_INLINE
3115 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3116             typename self_type::
3117 #endif
3118             const_reverse_iterator2 crend () const {
3119                 return rend ();
3120             }
3121 #endif
3122
3123             // Indices
3124             BOOST_UBLAS_INLINE
3125             size_type index1 () const {
3126                 return it2_.index1 ();
3127             }
3128             BOOST_UBLAS_INLINE
3129             size_type index2 () const {
3130                 return it2_.index2 ();
3131             }
3132
3133             // Assignment
3134             BOOST_UBLAS_INLINE
3135             const_iterator1 &operator = (const const_iterator1 &it) {
3136                 container_const_reference<self_type>::assign (&it ());
3137                 it1_ = it.it1_;
3138                 it2_ = it.it2_;
3139                 return *this;
3140             }
3141
3142             // Comparison
3143             BOOST_UBLAS_INLINE
3144             bool operator == (const const_iterator1 &it) const {
3145                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3146                 // FIXME we shouldn't compare floats
3147                 // BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
3148                 return it2_ == it.it2_;
3149             }
3150             BOOST_UBLAS_INLINE
3151             bool operator < (const const_iterator1 &it) const {
3152                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3153                 // FIXME we shouldn't compare floats
3154                 // BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
3155                 return it2_ < it.it2_;
3156             }
3157
3158         private:
3159             const_subiterator1_type it1_;
3160             const_iterator21_type it2_;
3161         };
3162 #endif
3163
3164         BOOST_UBLAS_INLINE
3165         const_iterator1 begin1 () const {
3166             return find1 (0, 0, 0);
3167         }
3168         BOOST_UBLAS_INLINE
3169         const_iterator1 cbegin1 () const {
3170             return begin1 ();
3171         }
3172         BOOST_UBLAS_INLINE
3173         const_iterator1 end1 () const {
3174             return find1 (0, size1 (), 0);
3175         }
3176         BOOST_UBLAS_INLINE
3177         const_iterator1 cend1 () const {
3178             return end1 ();
3179         }
3180
3181 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
3182         class const_iterator2:
3183             public container_const_reference<matrix_binary_scalar1>,
3184             public iterator_base_traits<typename E2::const_iterator2::iterator_category>::template
3185                 iterator_base<const_iterator2, value_type>::type {
3186         public:
3187             typedef typename E2::const_iterator2::iterator_category iterator_category;
3188             typedef typename matrix_binary_scalar1::difference_type difference_type;
3189             typedef typename matrix_binary_scalar1::value_type value_type;
3190             typedef typename matrix_binary_scalar1::const_reference reference;
3191             typedef typename matrix_binary_scalar1::const_pointer pointer;
3192
3193             typedef const_iterator1 dual_iterator_type;
3194             typedef const_reverse_iterator1 dual_reverse_iterator_type;
3195
3196             // Construction and destruction
3197             BOOST_UBLAS_INLINE
3198             const_iterator2 ():
3199                 container_const_reference<self_type> (), it1_ (), it2_ () {}
3200             BOOST_UBLAS_INLINE
3201             const_iterator2 (const self_type &mbs, const const_subiterator1_type &it1, const const_iterator22_type &it2):
3202                 container_const_reference<self_type> (mbs), it1_ (it1), it2_ (it2) {}
3203
3204             // Arithmetic
3205             BOOST_UBLAS_INLINE
3206             const_iterator2 &operator ++ () {
3207                 ++ it2_;
3208                 return *this;
3209             }
3210             BOOST_UBLAS_INLINE
3211             const_iterator2 &operator -- () {
3212                 -- it2_;
3213                 return *this;
3214             }
3215             BOOST_UBLAS_INLINE
3216             const_iterator2 &operator += (difference_type n) {
3217                 it2_ += n;
3218                 return *this;
3219             }
3220             BOOST_UBLAS_INLINE
3221             const_iterator2 &operator -= (difference_type n) {
3222                 it2_ -= n;
3223                 return *this;
3224             }
3225             BOOST_UBLAS_INLINE
3226             difference_type operator - (const const_iterator2 &it) const {
3227                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3228                 // FIXME we shouldn't compare floats
3229                 // BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
3230                 return it2_ - it.it2_;
3231             }
3232
3233             // Dereference
3234             BOOST_UBLAS_INLINE
3235             const_reference operator * () const {
3236                 return functor_type::apply (it1_, *it2_);
3237             }
3238             BOOST_UBLAS_INLINE
3239             const_reference operator [] (difference_type n) const {
3240                 return *(*this + n);
3241             }
3242
3243 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
3244             BOOST_UBLAS_INLINE
3245 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3246             typename self_type::
3247 #endif
3248             const_iterator1 begin () const {
3249                 return (*this) ().find1 (1, 0, index2 ());
3250             }
3251             BOOST_UBLAS_INLINE
3252 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3253             typename self_type::
3254 #endif
3255             const_iterator1 cbegin () const {
3256                 return begin ();
3257             }
3258             BOOST_UBLAS_INLINE
3259 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3260             typename self_type::
3261 #endif
3262             const_iterator1 end () const {
3263                 return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
3264             }
3265             BOOST_UBLAS_INLINE
3266 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3267             typename self_type::
3268 #endif
3269             const_iterator1 cend () const {
3270                 return end ();
3271             }
3272             BOOST_UBLAS_INLINE
3273 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3274             typename self_type::
3275 #endif
3276             const_reverse_iterator1 rbegin () const {
3277                 return const_reverse_iterator1 (end ());
3278             }
3279             BOOST_UBLAS_INLINE
3280 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3281             typename self_type::
3282 #endif
3283             const_reverse_iterator1 crbegin () const {
3284                 return rbegin ();
3285             }
3286             BOOST_UBLAS_INLINE
3287 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3288             typename self_type::
3289 #endif
3290             const_reverse_iterator1 rend () const {
3291                 return const_reverse_iterator1 (begin ());
3292             }
3293             BOOST_UBLAS_INLINE
3294 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3295             typename self_type::
3296 #endif
3297             const_reverse_iterator1 crend () const {
3298                 return rend ();
3299             }
3300 #endif
3301
3302             // Indices
3303             BOOST_UBLAS_INLINE
3304             size_type index1 () const {
3305                 return it2_.index1 ();
3306             }
3307             BOOST_UBLAS_INLINE
3308             size_type index2 () const {
3309                 return it2_.index2 ();
3310             }
3311
3312             // Assignment
3313             BOOST_UBLAS_INLINE
3314             const_iterator2 &operator = (const const_iterator2 &it) {
3315                 container_const_reference<self_type>::assign (&it ());
3316                 it1_ = it.it1_;
3317                 it2_ = it.it2_;
3318                 return *this;
3319             }
3320
3321             // Comparison
3322             BOOST_UBLAS_INLINE
3323             bool operator == (const const_iterator2 &it) const {
3324                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3325                 // FIXME we shouldn't compare floats
3326                 // BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
3327                 return it2_ == it.it2_;
3328             }
3329             BOOST_UBLAS_INLINE
3330             bool operator < (const const_iterator2 &it) const {
3331                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3332                 // FIXME we shouldn't compare floats
3333                 // BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
3334                 return it2_ < it.it2_;
3335             }
3336
3337         private:
3338             const_subiterator1_type it1_;
3339             const_iterator22_type it2_;
3340         };
3341 #endif
3342
3343         BOOST_UBLAS_INLINE
3344         const_iterator2 begin2 () const {
3345             return find2 (0, 0, 0);
3346         }
3347         BOOST_UBLAS_INLINE
3348         const_iterator2 cbegin2 () const {
3349             return begin2 ();
3350         }
3351         BOOST_UBLAS_INLINE
3352         const_iterator2 end2 () const {
3353             return find2 (0, 0, size2 ());
3354         }
3355         BOOST_UBLAS_INLINE
3356         const_iterator2 cend2 () const {
3357             return end2 ();
3358         }
3359
3360         // Reverse iterators
3361
3362         BOOST_UBLAS_INLINE
3363         const_reverse_iterator1 rbegin1 () const {
3364             return const_reverse_iterator1 (end1 ());
3365         }
3366         BOOST_UBLAS_INLINE
3367         const_reverse_iterator1 crbegin1 () const {
3368             return rbegin1 ();
3369         }
3370         BOOST_UBLAS_INLINE
3371         const_reverse_iterator1 rend1 () const {
3372             return const_reverse_iterator1 (begin1 ());
3373         }
3374         BOOST_UBLAS_INLINE
3375         const_reverse_iterator1 crend1 () const {
3376             return rend1 ();
3377         }
3378
3379         BOOST_UBLAS_INLINE
3380         const_reverse_iterator2 rbegin2 () const {
3381             return const_reverse_iterator2 (end2 ());
3382         }
3383         BOOST_UBLAS_INLINE
3384         const_reverse_iterator2 crbegin2 () const {
3385             return rbegin2 ();
3386         }
3387         BOOST_UBLAS_INLINE
3388         const_reverse_iterator2 rend2 () const {
3389             return const_reverse_iterator2 (begin2 ());
3390         }
3391         BOOST_UBLAS_INLINE
3392         const_reverse_iterator2 crend2 () const {
3393             return rend2 ();
3394         }
3395
3396     private:
3397         expression1_closure_type e1_;
3398         expression2_closure_type e2_;
3399     };
3400
3401     template<class E1, class E2, class F>
3402     struct matrix_binary_scalar1_traits {
3403         typedef matrix_binary_scalar1<E1, E2, F> expression_type;   // allow E1 to be builtin type
3404 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
3405         typedef expression_type result_type;
3406 #else
3407         typedef typename E2::matrix_temporary_type result_type;
3408 #endif
3409     };
3410
3411     // (t * m) [i] [j] = t * m [i] [j]
3412     template<class T1, class E2>
3413     BOOST_UBLAS_INLINE
3414     typename enable_if< is_convertible<T1, typename E2::value_type >,
3415     typename matrix_binary_scalar1_traits<const T1, E2, scalar_multiplies<T1, typename E2::value_type> >::result_type
3416     >::type
3417     operator * (const T1 &e1,
3418                 const matrix_expression<E2> &e2) {
3419         typedef typename matrix_binary_scalar1_traits<const T1, E2, scalar_multiplies<T1, typename E2::value_type> >::expression_type expression_type;
3420         return expression_type (e1, e2 ());
3421     }
3422
3423
3424     template<class E1, class E2, class F>
3425     class matrix_binary_scalar2:
3426         public matrix_expression<matrix_binary_scalar2<E1, E2, F> > {
3427
3428         typedef E1 expression1_type;
3429         typedef E2 expression2_type;
3430         typedef F functor_type;
3431     public:
3432         typedef typename E1::const_closure_type expression1_closure_type;
3433         typedef const E2& expression2_closure_type;
3434     private:
3435         typedef matrix_binary_scalar2<E1, E2, F> self_type;
3436     public:
3437 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
3438         using matrix_expression<self_type>::operator ();
3439 #endif
3440         typedef typename E1::size_type size_type;
3441         typedef typename E1::difference_type difference_type;
3442         typedef typename F::result_type value_type;
3443         typedef value_type const_reference;
3444         typedef const_reference reference;
3445
3446         typedef const self_type const_closure_type;
3447         typedef const_closure_type closure_type;
3448         typedef typename E1::orientation_category orientation_category;
3449         typedef unknown_storage_tag storage_category;
3450
3451         // Construction and destruction
3452         BOOST_UBLAS_INLINE
3453         matrix_binary_scalar2 (const expression1_type &e1, const expression2_type &e2):
3454             e1_ (e1), e2_ (e2) {}
3455
3456         // Accessors
3457         BOOST_UBLAS_INLINE
3458         size_type size1 () const {
3459             return e1_.size1 ();
3460         }
3461         BOOST_UBLAS_INLINE
3462         size_type size2 () const {
3463             return e1_.size2 ();
3464         }
3465
3466     public:
3467         // Element access
3468         BOOST_UBLAS_INLINE
3469         const_reference operator () (size_type i, size_type j) const {
3470             return functor_type::apply (e1_ (i, j), expression2_type (e2_));
3471         }
3472
3473         // Closure comparison
3474         BOOST_UBLAS_INLINE
3475         bool same_closure (const matrix_binary_scalar2 &mbs2) const {
3476             return (*this).e1_.same_closure (mbs2.e1_) &&
3477                    &e2_ == &(mbs2.e2_);
3478         }
3479
3480         // Iterator types
3481     private:
3482         typedef typename E1::const_iterator1 const_iterator11_type;
3483         typedef typename E1::const_iterator2 const_iterator12_type;
3484         typedef expression2_type const_subiterator2_type;
3485         typedef const value_type *const_pointer;
3486
3487     public:
3488 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
3489         typedef indexed_const_iterator1<const_closure_type, typename const_iterator11_type::iterator_category> const_iterator1;
3490         typedef const_iterator1 iterator1;
3491         typedef indexed_const_iterator2<const_closure_type, typename const_iterator12_type::iterator_category> const_iterator2;
3492         typedef const_iterator2 iterator2;
3493 #else
3494         class const_iterator1;
3495         typedef const_iterator1 iterator1;
3496         class const_iterator2;
3497         typedef const_iterator2 iterator2;
3498 #endif
3499         typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
3500         typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
3501
3502         // Element lookup
3503         BOOST_UBLAS_INLINE
3504         const_iterator1 find1 (int rank, size_type i, size_type j) const {
3505             const_iterator11_type it11 (e1_.find1 (rank, i, j));
3506 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
3507             return const_iterator1 (*this, it11.index1 (), it11.index2 ());
3508 #else
3509             return const_iterator1 (*this, it11, const_subiterator2_type (e2_));
3510 #endif
3511         }
3512         BOOST_UBLAS_INLINE
3513         const_iterator2 find2 (int rank, size_type i, size_type j) const {
3514             const_iterator12_type it12 (e1_.find2 (rank, i, j));
3515 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
3516             return const_iterator2 (*this, it12.index1 (), it12.index2 ());
3517 #else
3518             return const_iterator2 (*this, it12, const_subiterator2_type (e2_));
3519 #endif
3520         }
3521
3522         // Iterators enhance the iterators of the referenced expression
3523         // with the binary functor.
3524
3525 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
3526         class const_iterator1:
3527             public container_const_reference<matrix_binary_scalar2>,
3528             public iterator_base_traits<typename E1::const_iterator1::iterator_category>::template
3529                 iterator_base<const_iterator1, value_type>::type {
3530         public:
3531             typedef typename E1::const_iterator1::iterator_category iterator_category;
3532             typedef typename matrix_binary_scalar2::difference_type difference_type;
3533             typedef typename matrix_binary_scalar2::value_type value_type;
3534             typedef typename matrix_binary_scalar2::const_reference reference;
3535             typedef typename matrix_binary_scalar2::const_pointer pointer;
3536
3537             typedef const_iterator2 dual_iterator_type;
3538             typedef const_reverse_iterator2 dual_reverse_iterator_type;
3539
3540             // Construction and destruction
3541             BOOST_UBLAS_INLINE
3542             const_iterator1 ():
3543                 container_const_reference<self_type> (), it1_ (), it2_ () {}
3544             BOOST_UBLAS_INLINE
3545             const_iterator1 (const self_type &mbs, const const_iterator11_type &it1, const const_subiterator2_type &it2):
3546                 container_const_reference<self_type> (mbs), it1_ (it1), it2_ (it2) {}
3547
3548             // Arithmetic
3549             BOOST_UBLAS_INLINE
3550             const_iterator1 &operator ++ () {
3551                 ++ it1_;
3552                 return *this;
3553             }
3554             BOOST_UBLAS_INLINE
3555             const_iterator1 &operator -- () {
3556                 -- it1_ ;
3557                 return *this;
3558             }
3559             BOOST_UBLAS_INLINE
3560             const_iterator1 &operator += (difference_type n) {
3561                 it1_ += n;
3562                 return *this;
3563             }
3564             BOOST_UBLAS_INLINE
3565             const_iterator1 &operator -= (difference_type n) {
3566                 it1_ -= n;
3567                 return *this;
3568             }
3569             BOOST_UBLAS_INLINE
3570             difference_type operator - (const const_iterator1 &it) const {
3571                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3572                 // FIXME we shouldn't compare floats
3573                 // BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
3574                 return it1_ - it.it1_;
3575             }
3576
3577             // Dereference
3578             BOOST_UBLAS_INLINE
3579             const_reference operator * () const {
3580                 return functor_type::apply (*it1_, it2_);
3581             }
3582             BOOST_UBLAS_INLINE
3583             const_reference operator [] (difference_type n) const {
3584                 return *(*this + n);
3585             }
3586
3587 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
3588             BOOST_UBLAS_INLINE
3589 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3590             typename self_type::
3591 #endif
3592             const_iterator2 begin () const {
3593                 return (*this) ().find2 (1, index1 (), 0);
3594             }
3595             BOOST_UBLAS_INLINE
3596 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3597             typename self_type::
3598 #endif
3599             const_iterator2 cbegin () const {
3600                 return begin ();
3601             }
3602             BOOST_UBLAS_INLINE
3603 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3604             typename self_type::
3605 #endif
3606             const_iterator2 end () const {
3607                 return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
3608             }
3609             BOOST_UBLAS_INLINE
3610 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3611             typename self_type::
3612 #endif
3613             const_iterator2 cend () const {
3614                 return end ();
3615             }
3616             BOOST_UBLAS_INLINE
3617 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3618             typename self_type::
3619 #endif
3620             const_reverse_iterator2 rbegin () const {
3621                 return const_reverse_iterator2 (end ());
3622             }
3623             BOOST_UBLAS_INLINE
3624 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3625             typename self_type::
3626 #endif
3627             const_reverse_iterator2 crbegin () const {
3628                 return rbegin ();
3629             }
3630             BOOST_UBLAS_INLINE
3631 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3632             typename self_type::
3633 #endif
3634             const_reverse_iterator2 rend () const {
3635                 return const_reverse_iterator2 (begin ());
3636             }
3637             BOOST_UBLAS_INLINE
3638 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3639             typename self_type::
3640 #endif
3641             const_reverse_iterator2 crend () const {
3642                 return rend ();
3643             }
3644 #endif
3645
3646             // Indices
3647             BOOST_UBLAS_INLINE
3648             size_type index1 () const {
3649                 return it1_.index1 ();
3650             }
3651             BOOST_UBLAS_INLINE
3652             size_type index2 () const {
3653                 return it1_.index2 ();
3654             }
3655
3656             // Assignment
3657             BOOST_UBLAS_INLINE
3658             const_iterator1 &operator = (const const_iterator1 &it) {
3659                 container_const_reference<self_type>::assign (&it ());
3660                 it1_ = it.it1_;
3661                 it2_ = it.it2_;
3662                 return *this;
3663             }
3664
3665             // Comparison
3666             BOOST_UBLAS_INLINE
3667             bool operator == (const const_iterator1 &it) const {
3668                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3669                 // FIXME we shouldn't compare floats
3670                 // BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
3671                 return it1_ == it.it1_;
3672             }
3673             BOOST_UBLAS_INLINE
3674             bool operator < (const const_iterator1 &it) const {
3675                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3676                 // FIXME we shouldn't compare floats
3677                 // BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
3678                 return it1_ < it.it1_;
3679             }
3680
3681         private:
3682             const_iterator11_type it1_;
3683             const_subiterator2_type it2_;
3684         };
3685 #endif
3686
3687         BOOST_UBLAS_INLINE
3688         const_iterator1 begin1 () const {
3689             return find1 (0, 0, 0);
3690         }
3691         BOOST_UBLAS_INLINE
3692         const_iterator1 cbegin1 () const {
3693             return begin1 ();
3694         }
3695         BOOST_UBLAS_INLINE
3696         const_iterator1 end1 () const {
3697             return find1 (0, size1 (), 0);
3698         }
3699         BOOST_UBLAS_INLINE
3700         const_iterator1 cend1 () const {
3701             return end1 ();
3702         }
3703
3704 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
3705         class const_iterator2:
3706             public container_const_reference<matrix_binary_scalar2>,
3707             public iterator_base_traits<typename E1::const_iterator2::iterator_category>::template
3708                 iterator_base<const_iterator2, value_type>::type {
3709         public:
3710             typedef typename E1::const_iterator2::iterator_category iterator_category;
3711             typedef typename matrix_binary_scalar2::difference_type difference_type;
3712             typedef typename matrix_binary_scalar2::value_type value_type;
3713             typedef typename matrix_binary_scalar2::const_reference reference;
3714             typedef typename matrix_binary_scalar2::const_pointer pointer;
3715
3716             typedef const_iterator1 dual_iterator_type;
3717             typedef const_reverse_iterator1 dual_reverse_iterator_type;
3718
3719             // Construction and destruction
3720             BOOST_UBLAS_INLINE
3721             const_iterator2 ():
3722                 container_const_reference<self_type> (), it1_ (), it2_ () {}
3723             BOOST_UBLAS_INLINE
3724             const_iterator2 (const self_type &mbs, const const_iterator12_type &it1, const const_subiterator2_type &it2):
3725                 container_const_reference<self_type> (mbs), it1_ (it1), it2_ (it2) {}
3726
3727             // Arithmetic
3728             BOOST_UBLAS_INLINE
3729             const_iterator2 &operator ++ () {
3730                 ++ it1_;
3731                 return *this;
3732             }
3733             BOOST_UBLAS_INLINE
3734             const_iterator2 &operator -- () {
3735                 -- it1_;
3736                 return *this;
3737             }
3738             BOOST_UBLAS_INLINE
3739             const_iterator2 &operator += (difference_type n) {
3740                 it1_ += n;
3741                 return *this;
3742             }
3743             BOOST_UBLAS_INLINE
3744             const_iterator2 &operator -= (difference_type n) {
3745                 it1_ -= n;
3746                 return *this;
3747             }
3748             BOOST_UBLAS_INLINE
3749             difference_type operator - (const const_iterator2 &it) const {
3750                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3751                 // FIXME we shouldn't compare floats
3752                 // BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
3753                 return it1_ - it.it1_;
3754             }
3755
3756             // Dereference
3757             BOOST_UBLAS_INLINE
3758             const_reference operator * () const {
3759                 return functor_type::apply (*it1_, it2_);
3760             }
3761             BOOST_UBLAS_INLINE
3762             const_reference operator [] (difference_type n) const {
3763                 return *(*this + n);
3764             }
3765
3766 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
3767             BOOST_UBLAS_INLINE
3768 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3769             typename self_type::
3770 #endif
3771             const_iterator1 begin () const {
3772                 return (*this) ().find1 (1, 0, index2 ());
3773             }
3774             BOOST_UBLAS_INLINE
3775 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3776             typename self_type::
3777 #endif
3778             const_iterator1 cbegin () const {
3779                 return begin ();
3780             }
3781             BOOST_UBLAS_INLINE
3782 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3783             typename self_type::
3784 #endif
3785             const_iterator1 end () const {
3786                 return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
3787             }
3788             BOOST_UBLAS_INLINE
3789 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3790             typename self_type::
3791 #endif
3792             const_iterator1 cend () const {
3793                 return end ();
3794             }
3795             BOOST_UBLAS_INLINE
3796 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3797             typename self_type::
3798 #endif
3799             const_reverse_iterator1 rbegin () const {
3800                 return const_reverse_iterator1 (end ());
3801             }
3802             BOOST_UBLAS_INLINE
3803 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3804             typename self_type::
3805 #endif
3806             const_reverse_iterator1 crbegin () const {
3807                 return rbegin ();
3808             }
3809             BOOST_UBLAS_INLINE
3810 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3811             typename self_type::
3812 #endif
3813             const_reverse_iterator1 rend () const {
3814                 return const_reverse_iterator1 (begin ());
3815             }
3816             BOOST_UBLAS_INLINE
3817 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3818             typename self_type::
3819 #endif
3820             const_reverse_iterator1 crend () const {
3821                 return rend ();
3822             }
3823 #endif
3824
3825             // Indices
3826             BOOST_UBLAS_INLINE
3827             size_type index1 () const {
3828                 return it1_.index1 ();
3829             }
3830             BOOST_UBLAS_INLINE
3831             size_type index2 () const {
3832                 return it1_.index2 ();
3833             }
3834
3835             // Assignment
3836             BOOST_UBLAS_INLINE
3837             const_iterator2 &operator = (const const_iterator2 &it) {
3838                 container_const_reference<self_type>::assign (&it ());
3839                 it1_ = it.it1_;
3840                 it2_ = it.it2_;
3841                 return *this;
3842             }
3843
3844             // Comparison
3845             BOOST_UBLAS_INLINE
3846             bool operator == (const const_iterator2 &it) const {
3847                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3848                 // FIXME we shouldn't compare floats
3849                 // BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
3850                 return it1_ == it.it1_;
3851             }
3852             BOOST_UBLAS_INLINE
3853             bool operator < (const const_iterator2 &it) const {
3854                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
3855                 // FIXME we shouldn't compare floats
3856                 // BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
3857                 return it1_ < it.it1_;
3858             }
3859
3860         private:
3861             const_iterator12_type it1_;
3862             const_subiterator2_type it2_;
3863         };
3864 #endif
3865
3866         BOOST_UBLAS_INLINE
3867         const_iterator2 begin2 () const {
3868             return find2 (0, 0, 0);
3869         }
3870         BOOST_UBLAS_INLINE
3871         const_iterator2 cbegin2 () const {
3872             return begin2 ();
3873         }
3874         BOOST_UBLAS_INLINE
3875         const_iterator2 end2 () const {
3876             return find2 (0, 0, size2 ());
3877         }
3878         BOOST_UBLAS_INLINE
3879         const_iterator2 cend2 () const {
3880             return end2 ();
3881         }
3882
3883         // Reverse iterators
3884
3885         BOOST_UBLAS_INLINE
3886         const_reverse_iterator1 rbegin1 () const {
3887             return const_reverse_iterator1 (end1 ());
3888         }
3889         BOOST_UBLAS_INLINE
3890         const_reverse_iterator1 crbegin1 () const {
3891             return rbegin1 ();
3892         }
3893         BOOST_UBLAS_INLINE
3894         const_reverse_iterator1 rend1 () const {
3895             return const_reverse_iterator1 (begin1 ());
3896         }
3897         BOOST_UBLAS_INLINE
3898         const_reverse_iterator1 crend1 () const {
3899             return rend1 ();
3900         }
3901
3902         BOOST_UBLAS_INLINE
3903         const_reverse_iterator2 rbegin2 () const {
3904             return const_reverse_iterator2 (end2 ());
3905         }
3906         BOOST_UBLAS_INLINE
3907         const_reverse_iterator2 crbegin2 () const {
3908             return rbegin2 ();
3909         }
3910         BOOST_UBLAS_INLINE
3911         const_reverse_iterator2 rend2 () const {
3912             return const_reverse_iterator2 (begin2 ());
3913         }
3914         BOOST_UBLAS_INLINE
3915         const_reverse_iterator2 crend2 () const {
3916             return rend2 ();
3917         }
3918
3919     private:
3920         expression1_closure_type e1_;
3921         expression2_closure_type e2_;
3922     };
3923
3924     template<class E1, class E2, class F>
3925     struct matrix_binary_scalar2_traits {
3926         typedef matrix_binary_scalar2<E1, E2, F> expression_type;   // allow E2 to be builtin type
3927 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
3928         typedef expression_type result_type;
3929 #else
3930         typedef typename E1::matrix_temporary_type result_type;
3931 #endif
3932     };
3933
3934     // (m * t) [i] [j] = m [i] [j] * t
3935     template<class E1, class T2>
3936     BOOST_UBLAS_INLINE
3937     typename enable_if< is_convertible<T2, typename E1::value_type>,
3938     typename matrix_binary_scalar2_traits<E1, const T2, scalar_multiplies<typename E1::value_type, T2> >::result_type
3939     >::type
3940     operator * (const matrix_expression<E1> &e1,
3941                 const T2 &e2) {
3942         typedef typename matrix_binary_scalar2_traits<E1, const T2, scalar_multiplies<typename E1::value_type, T2> >::expression_type expression_type;
3943         return expression_type (e1 (), e2);
3944     }
3945
3946     // (m / t) [i] [j] = m [i] [j] / t
3947     template<class E1, class T2>
3948     BOOST_UBLAS_INLINE
3949     typename matrix_binary_scalar2_traits<E1, const T2, scalar_divides<typename E1::value_type, T2> >::result_type
3950     operator / (const matrix_expression<E1> &e1,
3951                 const T2 &e2) {
3952         typedef typename matrix_binary_scalar2_traits<E1, const T2, scalar_divides<typename E1::value_type, T2> >::expression_type expression_type;
3953         return expression_type (e1 (), e2);
3954     }
3955
3956
3957     template<class E1, class E2, class F>
3958     class matrix_vector_binary1:
3959         public vector_expression<matrix_vector_binary1<E1, E2, F> > {
3960
3961     public:
3962         typedef E1 expression1_type;
3963         typedef E2 expression2_type;
3964     private:
3965         typedef F functor_type;
3966     public:
3967         typedef typename E1::const_closure_type expression1_closure_type;
3968         typedef typename E2::const_closure_type expression2_closure_type;
3969     private:
3970         typedef matrix_vector_binary1<E1, E2, F> self_type;
3971     public:
3972 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
3973         using vector_expression<self_type>::operator ();
3974 #endif
3975         static const unsigned complexity = 1;
3976         typedef typename promote_traits<typename E1::size_type, typename E2::size_type>::promote_type size_type;
3977         typedef typename promote_traits<typename E1::difference_type, typename E2::difference_type>::promote_type difference_type;
3978         typedef typename F::result_type value_type;
3979         typedef value_type const_reference;
3980         typedef const_reference reference;
3981         typedef const self_type const_closure_type;
3982         typedef const_closure_type closure_type;
3983         typedef unknown_storage_tag storage_category;
3984
3985         // Construction and destruction
3986         BOOST_UBLAS_INLINE
3987         matrix_vector_binary1 (const expression1_type &e1, const expression2_type &e2):
3988             e1_ (e1), e2_ (e2) {}
3989
3990         // Accessors
3991         BOOST_UBLAS_INLINE
3992         size_type size () const {
3993             return e1_.size1 ();
3994         }
3995
3996     public:
3997         // Expression accessors
3998         BOOST_UBLAS_INLINE
3999         const expression1_closure_type &expression1 () const {
4000             return e1_;
4001         }
4002         BOOST_UBLAS_INLINE
4003         const expression2_closure_type &expression2 () const {
4004             return e2_;
4005         }
4006
4007     public:
4008         // Element access
4009         BOOST_UBLAS_INLINE
4010         const_reference operator () (size_type i) const {
4011             return functor_type::apply (e1_, e2_, i);
4012         }
4013
4014         // Closure comparison
4015         BOOST_UBLAS_INLINE
4016         bool same_closure (const matrix_vector_binary1 &mvb1) const {
4017             return (*this).expression1 ().same_closure (mvb1.expression1 ()) &&
4018                    (*this).expression2 ().same_closure (mvb1.expression2 ());
4019         }
4020
4021         // Iterator types
4022     private:
4023         typedef typename E1::const_iterator1 const_subiterator1_type;
4024         typedef typename E2::const_iterator const_subiterator2_type;
4025         typedef const value_type *const_pointer;
4026
4027     public:
4028 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
4029         typedef indexed_const_iterator<const_closure_type, typename const_subiterator1_type::iterator_category> const_iterator;
4030         typedef const_iterator iterator;
4031 #else
4032         class const_iterator;
4033         typedef const_iterator iterator;
4034 #endif
4035
4036         // Element lookup
4037         BOOST_UBLAS_INLINE
4038         const_iterator find (size_type i) const {
4039 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
4040             const_subiterator1_type it1 (e1_.find1 (0, i, 0));
4041             return const_iterator (*this, it1.index1 ());
4042 #else
4043             return const_iterator (*this, e1_.find1 (0, i, 0));
4044 #endif
4045         }
4046
4047
4048 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
4049         class const_iterator:
4050             public container_const_reference<matrix_vector_binary1>,
4051             public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
4052                                                                           typename E2::const_iterator::iterator_category>::iterator_category>::template
4053                 iterator_base<const_iterator, value_type>::type {
4054         public:
4055             typedef typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
4056                                                       typename E2::const_iterator::iterator_category>::iterator_category iterator_category;
4057             typedef typename matrix_vector_binary1::difference_type difference_type;
4058             typedef typename matrix_vector_binary1::value_type value_type;
4059             typedef typename matrix_vector_binary1::const_reference reference;
4060             typedef typename matrix_vector_binary1::const_pointer pointer;
4061
4062             // Construction and destruction
4063 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4064             BOOST_UBLAS_INLINE
4065             const_iterator ():
4066                 container_const_reference<self_type> (), it1_ (), e2_begin_ (), e2_end_ () {}
4067             BOOST_UBLAS_INLINE
4068             const_iterator (const self_type &mvb, const const_subiterator1_type &it1):
4069                 container_const_reference<self_type> (mvb), it1_ (it1), e2_begin_ (mvb.expression2 ().begin ()), e2_end_ (mvb.expression2 ().end ()) {}
4070 #else
4071             BOOST_UBLAS_INLINE
4072             const_iterator ():
4073                 container_const_reference<self_type> (), it1_ () {}
4074             BOOST_UBLAS_INLINE
4075             const_iterator (const self_type &mvb, const const_subiterator1_type &it1):
4076                 container_const_reference<self_type> (mvb), it1_ (it1) {}
4077 #endif
4078
4079         private:
4080             // Dense random access specialization
4081             BOOST_UBLAS_INLINE
4082             value_type dereference (dense_random_access_iterator_tag) const {
4083                 const self_type &mvb = (*this) ();
4084 #ifdef BOOST_UBLAS_USE_INDEXING
4085                 return mvb (index ());
4086 #elif BOOST_UBLAS_USE_ITERATING
4087                 difference_type size = BOOST_UBLAS_SAME (mvb.expression1 ().size2 (), mvb.expression2 ().size ());
4088 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4089                 return functor_type::apply (size, it1_.begin (), e2_begin_);
4090 #else
4091                 return functor_type::apply (size, it1_.begin (), mvb.expression2 ().begin ());
4092 #endif
4093 #else
4094                 difference_type size = BOOST_UBLAS_SAME (mvb.expression1 ().size2 (), mvb.expression2 ().size ());
4095                 if (size >= BOOST_UBLAS_ITERATOR_THRESHOLD)
4096 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4097                     return functor_type::apply (size, it1_.begin (), e2_begin_);
4098 #else
4099                     return functor_type::apply (size, it1_.begin (), mvb.expression2 ().begin ());
4100 #endif
4101                 else
4102                     return mvb (index ());
4103 #endif
4104             }
4105
4106             // Packed bidirectional specialization
4107             BOOST_UBLAS_INLINE
4108             value_type dereference (packed_random_access_iterator_tag) const {
4109 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4110                 return functor_type::apply (it1_.begin (), it1_.end (), e2_begin_, e2_end_);
4111 #else
4112                 const self_type &mvb = (*this) ();
4113 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
4114                 return functor_type::apply (it1_.begin (), it1_.end (),
4115                                         mvb.expression2 ().begin (), mvb.expression2 ().end ());
4116 #else
4117                 return functor_type::apply (boost::numeric::ublas::begin (it1_, iterator1_tag ()),
4118                                         boost::numeric::ublas::end (it1_, iterator1_tag ()),
4119                                         mvb.expression2 ().begin (), mvb.expression2 ().end ());
4120 #endif
4121 #endif
4122             }
4123
4124             // Sparse bidirectional specialization
4125             BOOST_UBLAS_INLINE
4126             value_type dereference (sparse_bidirectional_iterator_tag) const {
4127 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4128                 return functor_type::apply (it1_.begin (), it1_.end (), e2_begin_, e2_end_, sparse_bidirectional_iterator_tag ());
4129 #else
4130                 const self_type &mvb = (*this) ();
4131 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
4132                 return functor_type::apply (it1_.begin (), it1_.end (),
4133                                         mvb.expression2 ().begin (), mvb.expression2 ().end (), sparse_bidirectional_iterator_tag ());
4134 #else
4135                 return functor_type::apply (boost::numeric::ublas::begin (it1_, iterator1_tag ()),
4136                                         boost::numeric::ublas::end (it1_, iterator1_tag ()),
4137                                         mvb.expression2 ().begin (), mvb.expression2 ().end (), sparse_bidirectional_iterator_tag ());
4138 #endif
4139 #endif
4140             }
4141
4142         public:
4143             // Arithmetic
4144             BOOST_UBLAS_INLINE
4145             const_iterator &operator ++ () {
4146                 ++ it1_;
4147                 return *this;
4148             }
4149             BOOST_UBLAS_INLINE
4150             const_iterator &operator -- () {
4151                 -- it1_;
4152                 return *this;
4153             }
4154             BOOST_UBLAS_INLINE
4155             const_iterator &operator += (difference_type n) {
4156                 it1_ += n;
4157                 return *this;
4158             }
4159             BOOST_UBLAS_INLINE
4160             const_iterator &operator -= (difference_type n) {
4161                 it1_ -= n;
4162                 return *this;
4163             }
4164             BOOST_UBLAS_INLINE
4165             difference_type operator - (const const_iterator &it) const {
4166                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
4167                 return it1_ - it.it1_;
4168             }
4169
4170             // Dereference
4171             BOOST_UBLAS_INLINE
4172             const_reference operator * () const {
4173                 return dereference (iterator_category ());
4174             }
4175             BOOST_UBLAS_INLINE
4176             const_reference operator [] (difference_type n) const {
4177                 return *(*this + n);
4178             }
4179
4180             // Index
4181             BOOST_UBLAS_INLINE
4182             size_type index () const {
4183                 return it1_.index1 ();
4184             }
4185
4186             // Assignment
4187             BOOST_UBLAS_INLINE
4188             const_iterator &operator = (const const_iterator &it) {
4189                 container_const_reference<self_type>::assign (&it ());
4190                 it1_ = it.it1_;
4191 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4192                 e2_begin_ = it.e2_begin_;
4193                 e2_end_ = it.e2_end_;
4194 #endif
4195                 return *this;
4196             }
4197
4198             // Comparison
4199             BOOST_UBLAS_INLINE
4200             bool operator == (const const_iterator &it) const {
4201                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
4202                 return it1_ == it.it1_;
4203             }
4204             BOOST_UBLAS_INLINE
4205             bool operator < (const const_iterator &it) const {
4206                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
4207                 return it1_ < it.it1_;
4208             }
4209
4210         private:
4211             const_subiterator1_type it1_;
4212 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4213             // Mutable due to assignment
4214             /* const */ const_subiterator2_type e2_begin_;
4215             /* const */ const_subiterator2_type e2_end_;
4216 #endif
4217         };
4218 #endif
4219
4220         BOOST_UBLAS_INLINE
4221         const_iterator begin () const {
4222             return find (0);
4223         }
4224         BOOST_UBLAS_INLINE
4225         const_iterator cbegin () const {
4226             return begin ();
4227         }
4228         BOOST_UBLAS_INLINE
4229         const_iterator end () const {
4230             return find (size ());
4231         }
4232         BOOST_UBLAS_INLINE
4233         const_iterator cend () const {
4234             return end ();
4235         }
4236
4237         // Reverse iterator
4238         typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
4239
4240         BOOST_UBLAS_INLINE
4241         const_reverse_iterator rbegin () const {
4242             return const_reverse_iterator (end ());
4243         }
4244         BOOST_UBLAS_INLINE
4245         const_reverse_iterator crbegin () const {
4246             return rbegin ();
4247         }
4248         BOOST_UBLAS_INLINE
4249         const_reverse_iterator rend () const {
4250             return const_reverse_iterator (begin ());
4251         }
4252         BOOST_UBLAS_INLINE
4253         const_reverse_iterator crend () const {
4254             return rend ();
4255         }
4256
4257     private:
4258         expression1_closure_type e1_;
4259         expression2_closure_type e2_;
4260     };
4261
4262     template<class T1, class E1, class T2, class E2>
4263     struct matrix_vector_binary1_traits {
4264         typedef unknown_storage_tag storage_category;
4265         typedef row_major_tag orientation_category;
4266         typedef typename promote_traits<T1, T2>::promote_type promote_type;
4267         typedef matrix_vector_binary1<E1, E2, matrix_vector_prod1<E1, E2, promote_type> > expression_type;
4268 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
4269         typedef expression_type result_type;
4270 #else
4271         typedef typename E1::vector_temporary_type result_type;
4272 #endif
4273     };
4274
4275     template<class E1, class E2>
4276     BOOST_UBLAS_INLINE
4277     typename matrix_vector_binary1_traits<typename E1::value_type, E1,
4278                                           typename E2::value_type, E2>::result_type
4279     prod (const matrix_expression<E1> &e1,
4280           const vector_expression<E2> &e2,
4281           unknown_storage_tag,
4282           row_major_tag) {
4283         typedef typename matrix_vector_binary1_traits<typename E1::value_type, E1,
4284                                                       typename E2::value_type, E2>::expression_type expression_type;
4285         return expression_type (e1 (), e2 ());
4286     }
4287
4288     // Dispatcher
4289     template<class E1, class E2>
4290     BOOST_UBLAS_INLINE
4291     typename matrix_vector_binary1_traits<typename E1::value_type, E1,
4292                                           typename E2::value_type, E2>::result_type
4293     prod (const matrix_expression<E1> &e1,
4294           const vector_expression<E2> &e2) {
4295         BOOST_STATIC_ASSERT (E2::complexity == 0);
4296         typedef typename matrix_vector_binary1_traits<typename E1::value_type, E1,
4297                                                       typename E2::value_type, E2>::storage_category storage_category;
4298         typedef typename matrix_vector_binary1_traits<typename E1::value_type, E1,
4299                                                       typename E2::value_type, E2>::orientation_category orientation_category;
4300         return prod (e1, e2, storage_category (), orientation_category ());
4301     }
4302
4303     template<class E1, class E2>
4304     BOOST_UBLAS_INLINE
4305     typename matrix_vector_binary1_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4306                                           typename type_traits<typename E2::value_type>::precision_type, E2>::result_type
4307     prec_prod (const matrix_expression<E1> &e1,
4308                const vector_expression<E2> &e2,
4309                unknown_storage_tag,
4310                row_major_tag) {
4311         typedef typename matrix_vector_binary1_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4312                                                       typename type_traits<typename E2::value_type>::precision_type, E2>::expression_type expression_type;
4313         return expression_type (e1 (), e2 ());
4314     }
4315
4316     // Dispatcher
4317     template<class E1, class E2>
4318     BOOST_UBLAS_INLINE
4319     typename matrix_vector_binary1_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4320                                           typename type_traits<typename E2::value_type>::precision_type, E2>::result_type
4321     prec_prod (const matrix_expression<E1> &e1,
4322                const vector_expression<E2> &e2) {
4323         BOOST_STATIC_ASSERT (E2::complexity == 0);
4324         typedef typename matrix_vector_binary1_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4325                                                       typename type_traits<typename E2::value_type>::precision_type, E2>::storage_category storage_category;
4326         typedef typename matrix_vector_binary1_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4327                                                       typename type_traits<typename E2::value_type>::precision_type, E2>::orientation_category orientation_category;
4328         return prec_prod (e1, e2, storage_category (), orientation_category ());
4329     }
4330
4331     template<class V, class E1, class E2>
4332     BOOST_UBLAS_INLINE
4333     V &
4334     prod (const matrix_expression<E1> &e1,
4335           const vector_expression<E2> &e2,
4336           V &v) {
4337         return v.assign (prod (e1, e2));
4338     }
4339
4340     template<class V, class E1, class E2>
4341     BOOST_UBLAS_INLINE
4342     V &
4343     prec_prod (const matrix_expression<E1> &e1,
4344                const vector_expression<E2> &e2,
4345                V &v) {
4346         return v.assign (prec_prod (e1, e2));
4347     }
4348
4349     template<class V, class E1, class E2>
4350     BOOST_UBLAS_INLINE
4351     V
4352     prod (const matrix_expression<E1> &e1,
4353           const vector_expression<E2> &e2) {
4354         return V (prod (e1, e2));
4355     }
4356
4357     template<class V, class E1, class E2>
4358     BOOST_UBLAS_INLINE
4359     V
4360     prec_prod (const matrix_expression<E1> &e1,
4361                const vector_expression<E2> &e2) {
4362         return V (prec_prod (e1, e2));
4363     }
4364
4365     template<class E1, class E2, class F>
4366     class matrix_vector_binary2:
4367         public vector_expression<matrix_vector_binary2<E1, E2, F> > {
4368
4369         typedef E1 expression1_type;
4370         typedef E2 expression2_type;
4371         typedef F functor_type;
4372     public:
4373         typedef typename E1::const_closure_type expression1_closure_type;
4374         typedef typename E2::const_closure_type expression2_closure_type;
4375     private:
4376         typedef matrix_vector_binary2<E1, E2, F> self_type;
4377     public:
4378 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
4379         using vector_expression<self_type>::operator ();
4380 #endif
4381         static const unsigned complexity = 1;
4382         typedef typename promote_traits<typename E1::size_type, typename E2::size_type>::promote_type size_type;
4383         typedef typename promote_traits<typename E1::difference_type, typename E2::difference_type>::promote_type difference_type;
4384         typedef typename F::result_type value_type;
4385         typedef value_type const_reference;
4386         typedef const_reference reference;
4387         typedef const self_type const_closure_type;
4388         typedef const_closure_type closure_type;
4389         typedef unknown_storage_tag storage_category;
4390
4391         // Construction and destruction
4392         BOOST_UBLAS_INLINE
4393         matrix_vector_binary2 (const expression1_type &e1, const expression2_type &e2):
4394             e1_ (e1), e2_ (e2) {}
4395
4396         // Accessors
4397         BOOST_UBLAS_INLINE
4398         size_type size () const {
4399             return e2_.size2 ();
4400         }
4401
4402     public:
4403         // Expression accessors
4404         BOOST_UBLAS_INLINE
4405         const expression1_closure_type &expression1 () const {
4406             return e1_;
4407         }
4408         BOOST_UBLAS_INLINE
4409         const expression2_closure_type &expression2 () const {
4410             return e2_;
4411         }
4412     public:
4413
4414         // Element access
4415         BOOST_UBLAS_INLINE
4416         const_reference operator () (size_type j) const {
4417             return functor_type::apply (e1_, e2_, j);
4418         }
4419
4420         // Closure comparison
4421         BOOST_UBLAS_INLINE
4422         bool same_closure (const matrix_vector_binary2 &mvb2) const {
4423             return (*this).expression1 ().same_closure (mvb2.expression1 ()) &&
4424                    (*this).expression2 ().same_closure (mvb2.expression2 ());
4425         }
4426
4427         // Iterator types
4428     private:
4429         typedef typename E1::const_iterator const_subiterator1_type;
4430         typedef typename E2::const_iterator2 const_subiterator2_type;
4431         typedef const value_type *const_pointer;
4432
4433     public:
4434 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
4435         typedef indexed_const_iterator<const_closure_type, typename const_subiterator2_type::iterator_category> const_iterator;
4436         typedef const_iterator iterator;
4437 #else
4438         class const_iterator;
4439         typedef const_iterator iterator;
4440 #endif
4441
4442         // Element lookup
4443         BOOST_UBLAS_INLINE
4444         const_iterator find (size_type j) const {
4445 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
4446             const_subiterator2_type it2 (e2_.find2 (0, 0, j));
4447             return const_iterator (*this, it2.index2 ());
4448 #else
4449             return const_iterator (*this, e2_.find2 (0, 0, j));
4450 #endif
4451         }
4452
4453
4454 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
4455         class const_iterator:
4456             public container_const_reference<matrix_vector_binary2>,
4457             public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
4458                                                                           typename E2::const_iterator2::iterator_category>::iterator_category>::template
4459                 iterator_base<const_iterator, value_type>::type {
4460         public:
4461             typedef typename iterator_restrict_traits<typename E1::const_iterator::iterator_category,
4462                                                       typename E2::const_iterator2::iterator_category>::iterator_category iterator_category;
4463             typedef typename matrix_vector_binary2::difference_type difference_type;
4464             typedef typename matrix_vector_binary2::value_type value_type;
4465             typedef typename matrix_vector_binary2::const_reference reference;
4466             typedef typename matrix_vector_binary2::const_pointer pointer;
4467
4468             // Construction and destruction
4469 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4470             BOOST_UBLAS_INLINE
4471             const_iterator ():
4472                 container_const_reference<self_type> (), it2_ (), e1_begin_ (), e1_end_ () {}
4473             BOOST_UBLAS_INLINE
4474             const_iterator (const self_type &mvb, const const_subiterator2_type &it2):
4475                 container_const_reference<self_type> (mvb), it2_ (it2), e1_begin_ (mvb.expression1 ().begin ()), e1_end_ (mvb.expression1 ().end ()) {}
4476 #else
4477             BOOST_UBLAS_INLINE
4478             const_iterator ():
4479                 container_const_reference<self_type> (), it2_ () {}
4480             BOOST_UBLAS_INLINE
4481             const_iterator (const self_type &mvb, const const_subiterator2_type &it2):
4482                 container_const_reference<self_type> (mvb), it2_ (it2) {}
4483 #endif
4484
4485         private:
4486             // Dense random access specialization
4487             BOOST_UBLAS_INLINE
4488             value_type dereference (dense_random_access_iterator_tag) const {
4489                 const self_type &mvb = (*this) ();
4490 #ifdef BOOST_UBLAS_USE_INDEXING
4491                 return mvb (index ());
4492 #elif BOOST_UBLAS_USE_ITERATING
4493                 difference_type size = BOOST_UBLAS_SAME (mvb.expression2 ().size1 (), mvb.expression1 ().size ());
4494 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4495                 return functor_type::apply (size, e1_begin_, it2_.begin ());
4496 #else
4497                 return functor_type::apply (size, mvb.expression1 ().begin (), it2_.begin ());
4498 #endif
4499 #else
4500                 difference_type size = BOOST_UBLAS_SAME (mvb.expression2 ().size1 (), mvb.expression1 ().size ());
4501                 if (size >= BOOST_UBLAS_ITERATOR_THRESHOLD)
4502 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4503                     return functor_type::apply (size, e1_begin_, it2_.begin ());
4504 #else
4505                     return functor_type::apply (size, mvb.expression1 ().begin (), it2_.begin ());
4506 #endif
4507                 else
4508                     return mvb (index ());
4509 #endif
4510             }
4511
4512             // Packed bidirectional specialization
4513             BOOST_UBLAS_INLINE
4514             value_type dereference (packed_random_access_iterator_tag) const {
4515 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4516                 return functor_type::apply (e1_begin_, e1_end_, it2_.begin (), it2_.end ());
4517 #else
4518                 const self_type &mvb = (*this) ();
4519 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
4520                 return functor_type::apply (mvb.expression1 ().begin (), mvb.expression1 ().end (),
4521                                         it2_.begin (), it2_.end ());
4522 #else
4523                 return functor_type::apply (mvb.expression1 ().begin (), mvb.expression1 ().end (),
4524                                         boost::numeric::ublas::begin (it2_, iterator2_tag ()),
4525                                         boost::numeric::ublas::end (it2_, iterator2_tag ()));
4526 #endif
4527 #endif
4528             }
4529
4530             // Sparse bidirectional specialization
4531             BOOST_UBLAS_INLINE
4532             value_type dereference (sparse_bidirectional_iterator_tag) const {
4533 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4534                 return functor_type::apply (e1_begin_, e1_end_, it2_.begin (), it2_.end (), sparse_bidirectional_iterator_tag ());
4535 #else
4536                 const self_type &mvb = (*this) ();
4537 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
4538                 return functor_type::apply (mvb.expression1 ().begin (), mvb.expression1 ().end (),
4539                                         it2_.begin (), it2_.end (), sparse_bidirectional_iterator_tag ());
4540 #else
4541                 return functor_type::apply (mvb.expression1 ().begin (), mvb.expression1 ().end (),
4542                                         boost::numeric::ublas::begin (it2_, iterator2_tag ()),
4543                                         boost::numeric::ublas::end (it2_, iterator2_tag ()), sparse_bidirectional_iterator_tag ());
4544 #endif
4545 #endif
4546             }
4547
4548         public:
4549             // Arithmetic
4550             BOOST_UBLAS_INLINE
4551             const_iterator &operator ++ () {
4552                 ++ it2_;
4553                 return *this;
4554             }
4555             BOOST_UBLAS_INLINE
4556             const_iterator &operator -- () {
4557                 -- it2_;
4558                 return *this;
4559             }
4560             BOOST_UBLAS_INLINE
4561             const_iterator &operator += (difference_type n) {
4562                 it2_ += n;
4563                 return *this;
4564             }
4565             BOOST_UBLAS_INLINE
4566             const_iterator &operator -= (difference_type n) {
4567                 it2_ -= n;
4568                 return *this;
4569             }
4570             BOOST_UBLAS_INLINE
4571             difference_type operator - (const const_iterator &it) const {
4572                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
4573                 return it2_ - it.it2_;
4574             }
4575
4576             // Dereference
4577             BOOST_UBLAS_INLINE
4578             const_reference operator * () const {
4579                 return dereference (iterator_category ());
4580             }
4581             BOOST_UBLAS_INLINE
4582             const_reference operator [] (difference_type n) const {
4583                 return *(*this + n);
4584             }
4585
4586             // Index
4587             BOOST_UBLAS_INLINE
4588             size_type index () const {
4589                 return it2_.index2 ();
4590             }
4591
4592             // Assignment
4593             BOOST_UBLAS_INLINE
4594             const_iterator &operator = (const const_iterator &it) {
4595                 container_const_reference<self_type>::assign (&it ());
4596                 it2_ = it.it2_;
4597 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4598                 e1_begin_ = it.e1_begin_;
4599                 e1_end_ = it.e1_end_;
4600 #endif
4601                 return *this;
4602             }
4603
4604             // Comparison
4605             BOOST_UBLAS_INLINE
4606             bool operator == (const const_iterator &it) const {
4607                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
4608                 return it2_ == it.it2_;
4609             }
4610             BOOST_UBLAS_INLINE
4611             bool operator < (const const_iterator &it) const {
4612                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
4613                 return it2_ < it.it2_;
4614             }
4615
4616         private:
4617             const_subiterator2_type it2_;
4618 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4619             // Mutable due to assignment
4620             /* const */ const_subiterator1_type e1_begin_;
4621             /* const */ const_subiterator1_type e1_end_;
4622 #endif
4623         };
4624 #endif
4625
4626         BOOST_UBLAS_INLINE
4627         const_iterator begin () const {
4628             return find (0);
4629         }
4630         BOOST_UBLAS_INLINE
4631         const_iterator cbegin () const {
4632             return begin ();
4633         }
4634         BOOST_UBLAS_INLINE
4635         const_iterator end () const {
4636             return find (size ());
4637         }
4638         BOOST_UBLAS_INLINE
4639         const_iterator cend () const {
4640             return end ();
4641         }
4642
4643         // Reverse iterator
4644         typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
4645
4646         BOOST_UBLAS_INLINE
4647         const_reverse_iterator rbegin () const {
4648             return const_reverse_iterator (end ());
4649         }
4650         BOOST_UBLAS_INLINE
4651         const_reverse_iterator crbegin () const {
4652             return rbegin ();
4653         }
4654         BOOST_UBLAS_INLINE
4655         const_reverse_iterator rend () const {
4656             return const_reverse_iterator (begin ());
4657         }
4658         BOOST_UBLAS_INLINE
4659         const_reverse_iterator crend () const {
4660             return rend ();
4661         }
4662
4663     private:
4664         expression1_closure_type e1_;
4665         expression2_closure_type e2_;
4666     };
4667
4668     template<class T1, class E1, class T2, class E2>
4669     struct matrix_vector_binary2_traits {
4670         typedef unknown_storage_tag storage_category;
4671         typedef column_major_tag orientation_category;
4672         typedef typename promote_traits<T1, T2>::promote_type promote_type;
4673         typedef matrix_vector_binary2<E1, E2, matrix_vector_prod2<E1, E2, promote_type> > expression_type;
4674 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
4675         typedef expression_type result_type;
4676 #else
4677         typedef typename E2::vector_temporary_type result_type;
4678 #endif
4679     };
4680
4681     template<class E1, class E2>
4682     BOOST_UBLAS_INLINE
4683     typename matrix_vector_binary2_traits<typename E1::value_type, E1,
4684                                           typename E2::value_type, E2>::result_type
4685     prod (const vector_expression<E1> &e1,
4686           const matrix_expression<E2> &e2,
4687           unknown_storage_tag,
4688           column_major_tag) {
4689         typedef typename matrix_vector_binary2_traits<typename E1::value_type, E1,
4690                                                       typename E2::value_type, E2>::expression_type expression_type;
4691         return expression_type (e1 (), e2 ());
4692     }
4693
4694     // Dispatcher
4695     template<class E1, class E2>
4696     BOOST_UBLAS_INLINE
4697     typename matrix_vector_binary2_traits<typename E1::value_type, E1,
4698                                           typename E2::value_type, E2>::result_type
4699     prod (const vector_expression<E1> &e1,
4700           const matrix_expression<E2> &e2) {
4701         BOOST_STATIC_ASSERT (E1::complexity == 0);
4702         typedef typename matrix_vector_binary2_traits<typename E1::value_type, E1,
4703                                                       typename E2::value_type, E2>::storage_category storage_category;
4704         typedef typename matrix_vector_binary2_traits<typename E1::value_type, E1,
4705                                                       typename E2::value_type, E2>::orientation_category orientation_category;
4706         return prod (e1, e2, storage_category (), orientation_category ());
4707     }
4708
4709     template<class E1, class E2>
4710     BOOST_UBLAS_INLINE
4711     typename matrix_vector_binary2_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4712                                           typename type_traits<typename E2::value_type>::precision_type, E2>::result_type
4713     prec_prod (const vector_expression<E1> &e1,
4714                const matrix_expression<E2> &e2,
4715                unknown_storage_tag,
4716                column_major_tag) {
4717         typedef typename matrix_vector_binary2_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4718                                                       typename type_traits<typename E2::value_type>::precision_type, E2>::expression_type expression_type;
4719         return expression_type (e1 (), e2 ());
4720     }
4721
4722     // Dispatcher
4723     template<class E1, class E2>
4724     BOOST_UBLAS_INLINE
4725     typename matrix_vector_binary2_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4726                                           typename type_traits<typename E2::value_type>::precision_type, E2>::result_type
4727     prec_prod (const vector_expression<E1> &e1,
4728                const matrix_expression<E2> &e2) {
4729         BOOST_STATIC_ASSERT (E1::complexity == 0);
4730         typedef typename matrix_vector_binary2_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4731                                                       typename type_traits<typename E2::value_type>::precision_type, E2>::storage_category storage_category;
4732         typedef typename matrix_vector_binary2_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
4733                                                       typename type_traits<typename E2::value_type>::precision_type, E2>::orientation_category orientation_category;
4734         return prec_prod (e1, e2, storage_category (), orientation_category ());
4735     }
4736
4737     template<class V, class E1, class E2>
4738     BOOST_UBLAS_INLINE
4739     V &
4740     prod (const vector_expression<E1> &e1,
4741           const matrix_expression<E2> &e2,
4742           V &v) {
4743         return v.assign (prod (e1, e2));
4744     }
4745
4746     template<class V, class E1, class E2>
4747     BOOST_UBLAS_INLINE
4748     V &
4749     prec_prod (const vector_expression<E1> &e1,
4750                const matrix_expression<E2> &e2,
4751                V &v) {
4752         return v.assign (prec_prod (e1, e2));
4753     }
4754
4755     template<class V, class E1, class E2>
4756     BOOST_UBLAS_INLINE
4757     V
4758     prod (const vector_expression<E1> &e1,
4759           const matrix_expression<E2> &e2) {
4760         return V (prod (e1, e2));
4761     }
4762
4763     template<class V, class E1, class E2>
4764     BOOST_UBLAS_INLINE
4765     V
4766     prec_prod (const vector_expression<E1> &e1,
4767                const matrix_expression<E2> &e2) {
4768         return V (prec_prod (e1, e2));
4769     }
4770
4771     template<class E1, class E2, class F>
4772     class matrix_matrix_binary:
4773         public matrix_expression<matrix_matrix_binary<E1, E2, F> > {
4774
4775     public:
4776         typedef E1 expression1_type;
4777         typedef E2 expression2_type;
4778     private:
4779         typedef F functor_type;
4780     public:
4781         typedef typename E1::const_closure_type expression1_closure_type;
4782         typedef typename E2::const_closure_type expression2_closure_type;
4783     private:
4784         typedef matrix_matrix_binary<E1, E2, F> self_type;
4785     public:
4786 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
4787         using matrix_expression<self_type>::operator ();
4788 #endif
4789         static const unsigned complexity = 1;
4790         typedef typename promote_traits<typename E1::size_type, typename E2::size_type>::promote_type size_type;
4791         typedef typename promote_traits<typename E1::difference_type, typename E2::difference_type>::promote_type difference_type;
4792         typedef typename F::result_type value_type;
4793         typedef value_type const_reference;
4794         typedef const_reference reference;
4795         typedef const self_type const_closure_type;
4796         typedef const_closure_type closure_type;
4797         typedef unknown_orientation_tag orientation_category;
4798         typedef unknown_storage_tag storage_category;
4799
4800         // Construction and destruction
4801         BOOST_UBLAS_INLINE
4802         matrix_matrix_binary (const expression1_type &e1, const expression2_type &e2):
4803             e1_ (e1), e2_ (e2) {}
4804
4805         // Accessors
4806         BOOST_UBLAS_INLINE
4807         size_type size1 () const {
4808             return e1_.size1 ();
4809         }
4810         BOOST_UBLAS_INLINE
4811         size_type size2 () const {
4812             return e2_.size2 ();
4813         }
4814
4815     public:
4816         // Expression accessors
4817         BOOST_UBLAS_INLINE
4818         const expression1_closure_type &expression1 () const {
4819             return e1_;
4820         }
4821         BOOST_UBLAS_INLINE
4822         const expression2_closure_type &expression2 () const {
4823             return e2_;
4824         }
4825
4826     public:
4827         // Element access
4828         BOOST_UBLAS_INLINE
4829         const_reference operator () (size_type i, size_type j) const {
4830             return functor_type::apply (e1_, e2_, i, j);
4831         }
4832
4833         // Closure comparison
4834         BOOST_UBLAS_INLINE
4835         bool same_closure (const matrix_matrix_binary &mmb) const {
4836             return (*this).expression1 ().same_closure (mmb.expression1 ()) &&
4837                    (*this).expression2 ().same_closure (mmb.expression2 ());
4838         }
4839
4840         // Iterator types
4841     private:
4842         typedef typename E1::const_iterator1 const_iterator11_type;
4843         typedef typename E1::const_iterator2 const_iterator12_type;
4844         typedef typename E2::const_iterator1 const_iterator21_type;
4845         typedef typename E2::const_iterator2 const_iterator22_type;
4846         typedef const value_type *const_pointer;
4847
4848     public:
4849 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
4850         typedef typename iterator_restrict_traits<typename const_iterator11_type::iterator_category,
4851                                                   typename const_iterator22_type::iterator_category>::iterator_category iterator_category;
4852         typedef indexed_const_iterator1<const_closure_type, iterator_category> const_iterator1;
4853         typedef const_iterator1 iterator1;
4854         typedef indexed_const_iterator2<const_closure_type, iterator_category> const_iterator2;
4855         typedef const_iterator2 iterator2;
4856 #else
4857         class const_iterator1;
4858         typedef const_iterator1 iterator1;
4859         class const_iterator2;
4860         typedef const_iterator2 iterator2;
4861 #endif
4862         typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
4863         typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
4864
4865         // Element lookup
4866         BOOST_UBLAS_INLINE
4867         const_iterator1 find1 (int /* rank */, size_type i, size_type j) const {
4868             // FIXME sparse matrix tests fail!
4869             // const_iterator11_type it11 (e1_.find1 (rank, i, 0));
4870             const_iterator11_type it11 (e1_.find1 (0, i, 0));
4871 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
4872             return const_iterator1 (*this, it11.index1 (), j);
4873 #else
4874             // FIXME sparse matrix tests fail!
4875             // const_iterator22_type it22 (e2_.find2 (rank, 0, j));
4876             const_iterator22_type it22 (e2_.find2 (0, 0, j));
4877             return const_iterator1 (*this, it11, it22);
4878 #endif
4879         }
4880         BOOST_UBLAS_INLINE
4881         const_iterator2 find2 (int /* rank */, size_type i, size_type j) const {
4882             // FIXME sparse matrix tests fail!
4883             // const_iterator22_type it22 (e2_.find2 (rank, 0, j));
4884             const_iterator22_type it22 (e2_.find2 (0, 0, j));
4885 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
4886             return const_iterator2 (*this, i, it22.index2 ());
4887 #else
4888             // FIXME sparse matrix tests fail!
4889             // const_iterator11_type it11 (e1_.find1 (rank, i, 0));
4890             const_iterator11_type it11 (e1_.find1 (0, i, 0));
4891             return const_iterator2 (*this, it11, it22);
4892 #endif
4893         }
4894
4895
4896 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
4897         class const_iterator1:
4898             public container_const_reference<matrix_matrix_binary>,
4899             public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
4900                                                                           typename E2::const_iterator2::iterator_category>::iterator_category>::template
4901                 iterator_base<const_iterator1, value_type>::type {
4902         public:
4903             typedef typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
4904                                                       typename E2::const_iterator2::iterator_category>::iterator_category iterator_category;
4905             typedef typename matrix_matrix_binary::difference_type difference_type;
4906             typedef typename matrix_matrix_binary::value_type value_type;
4907             typedef typename matrix_matrix_binary::const_reference reference;
4908             typedef typename matrix_matrix_binary::const_pointer pointer;
4909
4910             typedef const_iterator2 dual_iterator_type;
4911             typedef const_reverse_iterator2 dual_reverse_iterator_type;
4912
4913             // Construction and destruction
4914 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4915             BOOST_UBLAS_INLINE
4916             const_iterator1 ():
4917                 container_const_reference<self_type> (), it1_ (), it2_ (), it2_begin_ (), it2_end_ () {}
4918             BOOST_UBLAS_INLINE
4919             const_iterator1 (const self_type &mmb, const const_iterator11_type &it1, const const_iterator22_type &it2):
4920                 container_const_reference<self_type> (mmb), it1_ (it1), it2_ (it2), it2_begin_ (it2.begin ()), it2_end_ (it2.end ()) {}
4921 #else
4922             BOOST_UBLAS_INLINE
4923             const_iterator1 ():
4924                 container_const_reference<self_type> (), it1_ (), it2_ () {}
4925             BOOST_UBLAS_INLINE
4926             const_iterator1 (const self_type &mmb, const const_iterator11_type &it1, const const_iterator22_type &it2):
4927                 container_const_reference<self_type> (mmb), it1_ (it1), it2_ (it2) {}
4928 #endif
4929
4930         private:
4931             // Random access specialization
4932             BOOST_UBLAS_INLINE
4933             value_type dereference (dense_random_access_iterator_tag) const {
4934                 const self_type &mmb = (*this) ();
4935 #ifdef BOOST_UBLAS_USE_INDEXING
4936                 return mmb (index1 (), index2 ());
4937 #elif BOOST_UBLAS_USE_ITERATING
4938                 difference_type size = BOOST_UBLAS_SAME (mmb.expression1 ().size2 (), mmb.expression2 ().size1 ());
4939 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4940                 return functor_type::apply (size, it1_.begin (), it2_begin_);
4941 #else
4942                 return functor_type::apply (size, it1_.begin (), it2_.begin ());
4943 #endif
4944 #else
4945                 difference_type size = BOOST_UBLAS_SAME (mmb.expression1 ().size2 (), mmb.expression2 ().size1 ());
4946                 if (size >= BOOST_UBLAS_ITERATOR_THRESHOLD)
4947 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4948                     return functor_type::apply (size, it1_.begin (), it2_begin_);
4949 #else
4950                     return functor_type::apply (size, it1_.begin (), it2_.begin ());
4951 #endif
4952                 else
4953                     return mmb (index1 (), index2 ());
4954 #endif
4955             }
4956
4957             // Packed bidirectional specialization
4958             BOOST_UBLAS_INLINE
4959             value_type dereference (packed_random_access_iterator_tag) const {
4960 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4961                 return functor_type::apply (it1_.begin (), it1_.end (),
4962                                         it2_begin_, it2_end_, packed_random_access_iterator_tag ());
4963 #else
4964 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
4965                 return functor_type::apply (it1_.begin (), it1_.end (),
4966                                         it2_.begin (), it2_.end (), packed_random_access_iterator_tag ());
4967 #else
4968                 return functor_type::apply (boost::numeric::ublas::begin (it1_, iterator1_tag ()),
4969                                         boost::numeric::ublas::end (it1_, iterator1_tag ()),
4970                                         boost::numeric::ublas::begin (it2_, iterator2_tag ()),
4971                                         boost::numeric::ublas::end (it2_, iterator2_tag ()), packed_random_access_iterator_tag ());
4972 #endif
4973 #endif
4974             }
4975
4976             // Sparse bidirectional specialization
4977             BOOST_UBLAS_INLINE
4978             value_type dereference (sparse_bidirectional_iterator_tag) const {
4979 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
4980                 return functor_type::apply (it1_.begin (), it1_.end (),
4981                                         it2_begin_, it2_end_, sparse_bidirectional_iterator_tag ());
4982 #else
4983 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
4984                 return functor_type::apply (it1_.begin (), it1_.end (),
4985                                         it2_.begin (), it2_.end (), sparse_bidirectional_iterator_tag ());
4986 #else
4987                 return functor_type::apply (boost::numeric::ublas::begin (it1_, iterator1_tag ()),
4988                                         boost::numeric::ublas::end (it1_, iterator1_tag ()),
4989                                         boost::numeric::ublas::begin (it2_, iterator2_tag ()),
4990                                         boost::numeric::ublas::end (it2_, iterator2_tag ()), sparse_bidirectional_iterator_tag ());
4991 #endif
4992 #endif
4993             }
4994
4995         public:
4996             // Arithmetic
4997             BOOST_UBLAS_INLINE
4998             const_iterator1 &operator ++ () {
4999                 ++ it1_;
5000                 return *this;
5001             }
5002             BOOST_UBLAS_INLINE
5003             const_iterator1 &operator -- () {
5004                 -- it1_;
5005                 return *this;
5006             }
5007             BOOST_UBLAS_INLINE
5008             const_iterator1 &operator += (difference_type n) {
5009                 it1_ += n;
5010                 return *this;
5011             }
5012             BOOST_UBLAS_INLINE
5013             const_iterator1 &operator -= (difference_type n) {
5014                 it1_ -= n;
5015                 return *this;
5016             }
5017             BOOST_UBLAS_INLINE
5018             difference_type operator - (const const_iterator1 &it) const {
5019                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
5020                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
5021                 return it1_ - it.it1_;
5022             }
5023
5024             // Dereference
5025             BOOST_UBLAS_INLINE
5026             const_reference operator * () const {
5027                 return dereference (iterator_category ());
5028             }
5029             BOOST_UBLAS_INLINE
5030             const_reference operator [] (difference_type n) const {
5031                 return *(*this + n);
5032             }
5033
5034 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
5035             BOOST_UBLAS_INLINE
5036 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5037             typename self_type::
5038 #endif
5039             const_iterator2 begin () const {
5040                 return (*this) ().find2 (1, index1 (), 0);
5041             }
5042             BOOST_UBLAS_INLINE
5043 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5044             typename self_type::
5045 #endif
5046             const_iterator2 cbegin () const {
5047                 return begin ();
5048             }
5049             BOOST_UBLAS_INLINE
5050 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5051             typename self_type::
5052 #endif
5053             const_iterator2 end () const {
5054                 return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
5055             }
5056             BOOST_UBLAS_INLINE
5057 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5058             typename self_type::
5059 #endif
5060             const_iterator2 cend () const {
5061                 return end ();
5062             }
5063             BOOST_UBLAS_INLINE
5064 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5065             typename self_type::
5066 #endif
5067             const_reverse_iterator2 rbegin () const {
5068                 return const_reverse_iterator2 (end ());
5069             }
5070             BOOST_UBLAS_INLINE
5071 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5072             typename self_type::
5073 #endif
5074             const_reverse_iterator2 crbegin () const {
5075                 return rbegin ();
5076             }
5077             BOOST_UBLAS_INLINE
5078 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5079             typename self_type::
5080 #endif
5081             const_reverse_iterator2 rend () const {
5082                 return const_reverse_iterator2 (begin ());
5083             }
5084             BOOST_UBLAS_INLINE
5085 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5086             typename self_type::
5087 #endif
5088             const_reverse_iterator2 crend () const {
5089                 return rend ();
5090             }
5091 #endif
5092
5093             // Indices
5094             BOOST_UBLAS_INLINE
5095             size_type index1 () const {
5096                 return it1_.index1 ();
5097             }
5098             BOOST_UBLAS_INLINE
5099             size_type index2 () const {
5100                 return it2_.index2 ();
5101             }
5102
5103             // Assignment
5104             BOOST_UBLAS_INLINE
5105             const_iterator1 &operator = (const const_iterator1 &it) {
5106                 container_const_reference<self_type>::assign (&it ());
5107                 it1_ = it.it1_;
5108                 it2_ = it.it2_;
5109 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
5110                 it2_begin_ = it.it2_begin_;
5111                 it2_end_ = it.it2_end_;
5112 #endif
5113                 return *this;
5114             }
5115
5116             // Comparison
5117             BOOST_UBLAS_INLINE
5118             bool operator == (const const_iterator1 &it) const {
5119                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
5120                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
5121                 return it1_ == it.it1_;
5122             }
5123             BOOST_UBLAS_INLINE
5124             bool operator < (const const_iterator1 &it) const {
5125                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
5126                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
5127                 return it1_ < it.it1_;
5128             }
5129
5130         private:
5131             const_iterator11_type it1_;
5132             // Mutable due to assignment
5133             /* const */ const_iterator22_type it2_;
5134 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
5135             /* const */ const_iterator21_type it2_begin_;
5136             /* const */ const_iterator21_type it2_end_;
5137 #endif
5138         };
5139 #endif
5140
5141         BOOST_UBLAS_INLINE
5142         const_iterator1 begin1 () const {
5143             return find1 (0, 0, 0);
5144         }
5145         BOOST_UBLAS_INLINE
5146         const_iterator1 cbegin1 () const {
5147             return begin1 ();
5148         }
5149         BOOST_UBLAS_INLINE
5150         const_iterator1 end1 () const {
5151             return find1 (0, size1 (), 0);
5152         }
5153         BOOST_UBLAS_INLINE
5154         const_iterator1 cend1 () const {
5155             return end1 ();
5156         }
5157
5158 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
5159         class const_iterator2:
5160             public container_const_reference<matrix_matrix_binary>,
5161             public iterator_base_traits<typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
5162                                                                           typename E2::const_iterator2::iterator_category>::iterator_category>::template
5163                 iterator_base<const_iterator2, value_type>::type {
5164         public:
5165             typedef typename iterator_restrict_traits<typename E1::const_iterator1::iterator_category,
5166                                                       typename E2::const_iterator2::iterator_category>::iterator_category iterator_category;
5167             typedef typename matrix_matrix_binary::difference_type difference_type;
5168             typedef typename matrix_matrix_binary::value_type value_type;
5169             typedef typename matrix_matrix_binary::const_reference reference;
5170             typedef typename matrix_matrix_binary::const_pointer pointer;
5171
5172             typedef const_iterator1 dual_iterator_type;
5173             typedef const_reverse_iterator1 dual_reverse_iterator_type;
5174
5175             // Construction and destruction
5176 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
5177             BOOST_UBLAS_INLINE
5178             const_iterator2 ():
5179                 container_const_reference<self_type> (), it1_ (), it2_ (), it1_begin_ (), it1_end_ () {}
5180             BOOST_UBLAS_INLINE
5181             const_iterator2 (const self_type &mmb, const const_iterator11_type &it1, const const_iterator22_type &it2):
5182                 container_const_reference<self_type> (mmb), it1_ (it1), it2_ (it2), it1_begin_ (it1.begin ()), it1_end_ (it1.end ()) {}
5183 #else
5184             BOOST_UBLAS_INLINE
5185             const_iterator2 ():
5186                 container_const_reference<self_type> (), it1_ (), it2_ () {}
5187             BOOST_UBLAS_INLINE
5188             const_iterator2 (const self_type &mmb, const const_iterator11_type &it1, const const_iterator22_type &it2):
5189                 container_const_reference<self_type> (mmb), it1_ (it1), it2_ (it2) {}
5190 #endif
5191
5192         private:
5193             // Random access specialization
5194             BOOST_UBLAS_INLINE
5195             value_type dereference (dense_random_access_iterator_tag) const {
5196                 const self_type &mmb = (*this) ();
5197 #ifdef BOOST_UBLAS_USE_INDEXING
5198                 return mmb (index1 (), index2 ());
5199 #elif BOOST_UBLAS_USE_ITERATING
5200                 difference_type size = BOOST_UBLAS_SAME (mmb.expression1 ().size2 (), mmb.expression2 ().size1 ());
5201 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
5202                 return functor_type::apply (size, it1_begin_, it2_.begin ());
5203 #else
5204                 return functor_type::apply (size, it1_.begin (), it2_.begin ());
5205 #endif
5206 #else
5207                 difference_type size = BOOST_UBLAS_SAME (mmb.expression1 ().size2 (), mmb.expression2 ().size1 ());
5208                 if (size >= BOOST_UBLAS_ITERATOR_THRESHOLD)
5209 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
5210                     return functor_type::apply (size, it1_begin_, it2_.begin ());
5211 #else
5212                     return functor_type::apply (size, it1_.begin (), it2_.begin ());
5213 #endif
5214                 else
5215                     return mmb (index1 (), index2 ());
5216 #endif
5217             }
5218
5219             // Packed bidirectional specialization
5220             BOOST_UBLAS_INLINE
5221             value_type dereference (packed_random_access_iterator_tag) const {
5222 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
5223                 return functor_type::apply (it1_begin_, it1_end_,
5224                                         it2_.begin (), it2_.end (), packed_random_access_iterator_tag ());
5225 #else
5226 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
5227                 return functor_type::apply (it1_.begin (), it1_.end (),
5228                                         it2_.begin (), it2_.end (), packed_random_access_iterator_tag ());
5229 #else
5230                 return functor_type::apply (boost::numeric::ublas::begin (it1_, iterator1_tag ()),
5231                                         boost::numeric::ublas::end (it1_, iterator1_tag ()),
5232                                         boost::numeric::ublas::begin (it2_, iterator2_tag ()),
5233                                         boost::numeric::ublas::end (it2_, iterator2_tag ()), packed_random_access_iterator_tag ());
5234 #endif
5235 #endif
5236             }
5237
5238             // Sparse bidirectional specialization
5239             BOOST_UBLAS_INLINE
5240             value_type dereference (sparse_bidirectional_iterator_tag) const {
5241 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
5242                 return functor_type::apply (it1_begin_, it1_end_,
5243                                         it2_.begin (), it2_.end (), sparse_bidirectional_iterator_tag ());
5244 #else
5245 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
5246                 return functor_type::apply (it1_.begin (), it1_.end (),
5247                                         it2_.begin (), it2_.end (), sparse_bidirectional_iterator_tag ());
5248 #else
5249                 return functor_type::apply (boost::numeric::ublas::begin (it1_, iterator1_tag ()),
5250                                         boost::numeric::ublas::end (it1_, iterator1_tag ()),
5251                                         boost::numeric::ublas::begin (it2_, iterator2_tag ()),
5252                                         boost::numeric::ublas::end (it2_, iterator2_tag ()), sparse_bidirectional_iterator_tag ());
5253 #endif
5254 #endif
5255             }
5256
5257         public:
5258             // Arithmetic
5259             BOOST_UBLAS_INLINE
5260             const_iterator2 &operator ++ () {
5261                 ++ it2_;
5262                 return *this;
5263             }
5264             BOOST_UBLAS_INLINE
5265             const_iterator2 &operator -- () {
5266                 -- it2_;
5267                 return *this;
5268             }
5269             BOOST_UBLAS_INLINE
5270             const_iterator2 &operator += (difference_type n) {
5271                 it2_ += n;
5272                 return *this;
5273             }
5274             BOOST_UBLAS_INLINE
5275             const_iterator2 &operator -= (difference_type n) {
5276                 it2_ -= n;
5277                 return *this;
5278             }
5279             BOOST_UBLAS_INLINE
5280             difference_type operator - (const const_iterator2 &it) const {
5281                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
5282                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
5283                 return it2_ - it.it2_;
5284             }
5285
5286             // Dereference
5287             BOOST_UBLAS_INLINE
5288             const_reference operator * () const {
5289                 return dereference (iterator_category ());
5290             }
5291             BOOST_UBLAS_INLINE
5292             const_reference operator [] (difference_type n) const {
5293                 return *(*this + n);
5294             }
5295
5296 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
5297             BOOST_UBLAS_INLINE
5298 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5299             typename self_type::
5300 #endif
5301             const_iterator1 begin () const {
5302                 return (*this) ().find1 (1, 0, index2 ());
5303             }
5304             BOOST_UBLAS_INLINE
5305 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5306             typename self_type::
5307 #endif
5308             const_iterator1 cbegin () const {
5309                 return begin ();
5310             }
5311             BOOST_UBLAS_INLINE
5312 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5313             typename self_type::
5314 #endif
5315             const_iterator1 end () const {
5316                 return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
5317             }
5318             BOOST_UBLAS_INLINE
5319 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5320             typename self_type::
5321 #endif
5322             const_iterator1 cend () const {
5323                 return end ();
5324             }
5325             BOOST_UBLAS_INLINE
5326 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5327             typename self_type::
5328 #endif
5329             const_reverse_iterator1 rbegin () const {
5330                 return const_reverse_iterator1 (end ());
5331             }
5332             BOOST_UBLAS_INLINE
5333 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5334             typename self_type::
5335 #endif
5336             const_reverse_iterator1 crbegin () const {
5337                 return rbegin ();
5338             }
5339             BOOST_UBLAS_INLINE
5340 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5341             typename self_type::
5342 #endif
5343             const_reverse_iterator1 rend () const {
5344                 return const_reverse_iterator1 (begin ());
5345             }
5346             BOOST_UBLAS_INLINE
5347 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5348             typename self_type::
5349 #endif
5350             const_reverse_iterator1 crend () const {
5351                 return rend ();
5352             }
5353 #endif
5354
5355             // Indices
5356             BOOST_UBLAS_INLINE
5357             size_type index1 () const {
5358                 return it1_.index1 ();
5359             }
5360             BOOST_UBLAS_INLINE
5361             size_type index2 () const {
5362                 return it2_.index2 ();
5363             }
5364
5365             // Assignment
5366             BOOST_UBLAS_INLINE
5367             const_iterator2 &operator = (const const_iterator2 &it) {
5368                 container_const_reference<self_type>::assign (&it ());
5369                 it1_ = it.it1_;
5370                 it2_ = it.it2_;
5371 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
5372                 it1_begin_ = it.it1_begin_;
5373                 it1_end_ = it.it1_end_;
5374 #endif
5375                 return *this;
5376             }
5377
5378             // Comparison
5379             BOOST_UBLAS_INLINE
5380             bool operator == (const const_iterator2 &it) const {
5381                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
5382                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
5383                 return it2_ == it.it2_;
5384             }
5385             BOOST_UBLAS_INLINE
5386             bool operator < (const const_iterator2 &it) const {
5387                 BOOST_UBLAS_CHECK ((*this) ().same_closure (it ()), external_logic ());
5388                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
5389                 return it2_ < it.it2_;
5390             }
5391
5392         private:
5393             // Mutable due to assignment
5394             /* const */ const_iterator11_type it1_;
5395             const_iterator22_type it2_;
5396 #ifdef BOOST_UBLAS_USE_INVARIANT_HOISTING
5397             /* const */ const_iterator12_type it1_begin_;
5398             /* const */ const_iterator12_type it1_end_;
5399 #endif
5400         };
5401 #endif
5402
5403         BOOST_UBLAS_INLINE
5404         const_iterator2 begin2 () const {
5405             return find2 (0, 0, 0);
5406         }
5407         BOOST_UBLAS_INLINE
5408         const_iterator2 cbegin2 () const {
5409             return begin2 ();
5410         }
5411         BOOST_UBLAS_INLINE
5412         const_iterator2 end2 () const {
5413             return find2 (0, 0, size2 ());
5414         }
5415         BOOST_UBLAS_INLINE
5416         const_iterator2 cend2 () const {
5417             return end2 ();
5418         }
5419
5420         // Reverse iterators
5421
5422         BOOST_UBLAS_INLINE
5423         const_reverse_iterator1 rbegin1 () const {
5424             return const_reverse_iterator1 (end1 ());
5425         }
5426         BOOST_UBLAS_INLINE
5427         const_reverse_iterator1 crbegin1 () const {
5428             return rbegin1 ();
5429         }
5430         BOOST_UBLAS_INLINE
5431         const_reverse_iterator1 rend1 () const {
5432             return const_reverse_iterator1 (begin1 ());
5433         }
5434         BOOST_UBLAS_INLINE
5435         const_reverse_iterator1 crend1 () const {
5436             return rend1 ();
5437         }
5438
5439         BOOST_UBLAS_INLINE
5440         const_reverse_iterator2 rbegin2 () const {
5441             return const_reverse_iterator2 (end2 ());
5442         }
5443         BOOST_UBLAS_INLINE
5444         const_reverse_iterator2 crbegin2 () const {
5445             return rbegin2 ();
5446         }
5447         BOOST_UBLAS_INLINE
5448         const_reverse_iterator2 rend2 () const {
5449             return const_reverse_iterator2 (begin2 ());
5450         }
5451         BOOST_UBLAS_INLINE
5452         const_reverse_iterator2 crend2 () const {
5453             return rend2 ();
5454         }
5455
5456     private:
5457         expression1_closure_type e1_;
5458         expression2_closure_type e2_;
5459     };
5460
5461     template<class T1, class E1, class T2, class E2>
5462     struct matrix_matrix_binary_traits {
5463         typedef unknown_storage_tag storage_category;
5464         typedef unknown_orientation_tag orientation_category;
5465         typedef typename promote_traits<T1, T2>::promote_type promote_type;
5466         typedef matrix_matrix_binary<E1, E2, matrix_matrix_prod<E1, E2, promote_type> > expression_type;
5467 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
5468         typedef expression_type result_type;
5469 #else
5470         typedef typename E1::matrix_temporary_type result_type;
5471 #endif
5472     };
5473
5474     template<class E1, class E2>
5475     BOOST_UBLAS_INLINE
5476     typename matrix_matrix_binary_traits<typename E1::value_type, E1,
5477                                          typename E2::value_type, E2>::result_type
5478     prod (const matrix_expression<E1> &e1,
5479           const matrix_expression<E2> &e2,
5480           unknown_storage_tag,
5481           unknown_orientation_tag) {
5482         typedef typename matrix_matrix_binary_traits<typename E1::value_type, E1,
5483                                                      typename E2::value_type, E2>::expression_type expression_type;
5484         return expression_type (e1 (), e2 ());
5485     }
5486
5487     // Dispatcher
5488     template<class E1, class E2>
5489     BOOST_UBLAS_INLINE
5490     typename matrix_matrix_binary_traits<typename E1::value_type, E1,
5491                                          typename E2::value_type, E2>::result_type
5492     prod (const matrix_expression<E1> &e1,
5493           const matrix_expression<E2> &e2) {
5494         BOOST_STATIC_ASSERT (E1::complexity == 0 && E2::complexity == 0);
5495         typedef typename matrix_matrix_binary_traits<typename E1::value_type, E1,
5496                                                      typename E2::value_type, E2>::storage_category storage_category;
5497         typedef typename matrix_matrix_binary_traits<typename E1::value_type, E1,
5498                                                      typename E2::value_type, E2>::orientation_category orientation_category;
5499         return prod (e1, e2, storage_category (), orientation_category ());
5500     }
5501
5502     template<class E1, class E2>
5503     BOOST_UBLAS_INLINE
5504     typename matrix_matrix_binary_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
5505                                          typename type_traits<typename E2::value_type>::precision_type, E2>::result_type
5506     prec_prod (const matrix_expression<E1> &e1,
5507                const matrix_expression<E2> &e2,
5508                unknown_storage_tag,
5509                unknown_orientation_tag) {
5510         typedef typename matrix_matrix_binary_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
5511                                                      typename type_traits<typename E2::value_type>::precision_type, E2>::expression_type expression_type;
5512         return expression_type (e1 (), e2 ());
5513     }
5514
5515     // Dispatcher
5516     template<class E1, class E2>
5517     BOOST_UBLAS_INLINE
5518     typename matrix_matrix_binary_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
5519                                          typename type_traits<typename E2::value_type>::precision_type, E2>::result_type
5520     prec_prod (const matrix_expression<E1> &e1,
5521                const matrix_expression<E2> &e2) {
5522         BOOST_STATIC_ASSERT (E1::complexity == 0 && E2::complexity == 0);
5523         typedef typename matrix_matrix_binary_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
5524                                                      typename type_traits<typename E2::value_type>::precision_type, E2>::storage_category storage_category;
5525         typedef typename matrix_matrix_binary_traits<typename type_traits<typename E1::value_type>::precision_type, E1,
5526                                                      typename type_traits<typename E2::value_type>::precision_type, E2>::orientation_category orientation_category;
5527         return prec_prod (e1, e2, storage_category (), orientation_category ());
5528     }
5529
5530     template<class M, class E1, class E2>
5531     BOOST_UBLAS_INLINE
5532     M &
5533     prod (const matrix_expression<E1> &e1,
5534           const matrix_expression<E2> &e2,
5535           M &m) {
5536         return m.assign (prod (e1, e2));
5537     }
5538
5539     template<class M, class E1, class E2>
5540     BOOST_UBLAS_INLINE
5541     M &
5542     prec_prod (const matrix_expression<E1> &e1,
5543                const matrix_expression<E2> &e2,
5544                M &m) {
5545         return m.assign (prec_prod (e1, e2));
5546     }
5547
5548     template<class M, class E1, class E2>
5549     BOOST_UBLAS_INLINE
5550     M
5551     prod (const matrix_expression<E1> &e1,
5552           const matrix_expression<E2> &e2) {
5553         return M (prod (e1, e2));
5554     }
5555
5556     template<class M, class E1, class E2>
5557     BOOST_UBLAS_INLINE
5558     M
5559     prec_prod (const matrix_expression<E1> &e1,
5560                const matrix_expression<E2> &e2) {
5561         return M (prec_prod (e1, e2));
5562     }
5563
5564     template<class E, class F>
5565     class matrix_scalar_unary:
5566         public scalar_expression<matrix_scalar_unary<E, F> > {
5567     public:
5568         typedef E expression_type;
5569         typedef F functor_type;
5570         typedef typename F::result_type value_type;
5571         typedef typename E::const_closure_type expression_closure_type;
5572
5573         // Construction and destruction
5574         BOOST_UBLAS_INLINE
5575         explicit matrix_scalar_unary (const expression_type &e):
5576             e_ (e) {}
5577
5578     private:
5579         // Expression accessors
5580         BOOST_UBLAS_INLINE
5581         const expression_closure_type &expression () const {
5582             return e_;
5583         }
5584
5585     public:
5586         BOOST_UBLAS_INLINE
5587         operator value_type () const {
5588             return functor_type::apply (e_);
5589         }
5590
5591     private:
5592         expression_closure_type e_;
5593     };
5594
5595     template<class E, class F>
5596     struct matrix_scalar_unary_traits {
5597         typedef matrix_scalar_unary<E, F> expression_type;
5598 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
5599          typedef expression_type result_type;
5600 #else
5601          typedef typename F::result_type result_type;
5602 #endif
5603     };
5604
5605     template<class E>
5606     BOOST_UBLAS_INLINE
5607     typename matrix_scalar_unary_traits<E, matrix_norm_1<E> >::result_type
5608     norm_1 (const matrix_expression<E> &e) {
5609         typedef typename matrix_scalar_unary_traits<E, matrix_norm_1<E> >::expression_type expression_type;
5610         return expression_type (e ());
5611     }
5612
5613     template<class E>
5614     BOOST_UBLAS_INLINE
5615     typename matrix_scalar_unary_traits<E, matrix_norm_frobenius<E> >::result_type
5616     norm_frobenius (const matrix_expression<E> &e) {
5617         typedef typename matrix_scalar_unary_traits<E, matrix_norm_frobenius<E> >::expression_type expression_type;
5618         return expression_type (e ());
5619     }
5620
5621     template<class E>
5622     BOOST_UBLAS_INLINE
5623     typename matrix_scalar_unary_traits<E, matrix_norm_inf<E> >::result_type
5624     norm_inf (const matrix_expression<E> &e) {
5625         typedef typename matrix_scalar_unary_traits<E, matrix_norm_inf<E> >::expression_type expression_type;
5626         return expression_type (e ());
5627     }
5628
5629 }}}
5630
5631 #endif