ceph/src/boost/libs/numeric/ublas/include/boost/numeric/ublas/matrix.hpp

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