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

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