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

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