boost/numeric/ublas/matrix_proxy.hpp

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