source: NonGTP/Boost/boost/numeric/ublas/matrix_expression.hpp @ 857

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