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

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