New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
matrix.hpp in vendors/XIOS/current/extern/boost/include/boost/numeric/ublas – NEMO

source: vendors/XIOS/current/extern/boost/include/boost/numeric/ublas/matrix.hpp @ 3428

Last change on this file since 3428 was 3428, checked in by rblod, 12 years ago

importing initial XIOS vendor drop

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