Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:43:11

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