source: vendor/nemo/current/NEMOGCM/EXTERNAL/XIOS/extern/boost/include/boost/numeric/ublas/matrix_sparse.hpp @ 44

Last change on this file since 44 was 44, checked in by cholod, 12 years ago

Load NEMO_TMP into vendor/nemo/current.

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