- Main Page
- Related Pages
- Classes
- Files

- File List

# matrix_sparse.hpp

    00001 //
    00002 //  Copyright (c) 2000-2007
    00003 //  Joerg Walter, Mathias Koch, Gunter Winkler
    00004 //
    00005 //  Distributed under the Boost Software License, Version 1.0. (See
    00006 //  accompanying file LICENSE_1_0.txt or copy at
    00007 //  http://www.boost.org/LICENSE_1_0.txt)
    00008 //
    00009 //  The authors gratefully acknowledge the support of
    00010 //  GeNeSys mbH & Co. KG in producing this work.
    00011 //
    00012
    00013 #ifndef _BOOST_UBLAS_MATRIX_SPARSE_
    00014 #define _BOOST_UBLAS_MATRIX_SPARSE_
    00015
    00016 #include <boost/numeric/ublas/vector_sparse.hpp>
    00017 #include <boost/numeric/ublas/matrix_expression.hpp>
    00018 #include <boost/numeric/ublas/detail/matrix_assign.hpp>
    00019 #if BOOST_UBLAS_TYPE_CHECK
    00020 #include <boost/numeric/ublas/matrix.hpp>
    00021 #endif
    00022
    00023 // Iterators based on ideas of Jeremy Siek
    00024
    00025 namespace boost { namespace numeric { namespace ublas {
    00026
    00027 #ifdef BOOST_UBLAS_STRICT_MATRIX_SPARSE
    00028
    00029     template<class M>
    00030     class sparse_matrix_element:
    00031        public container_reference<M> {
    00032     public:
    00033         typedef M matrix_type;
    00034         typedef typename M::size_type size_type;
    00035         typedef typename M::value_type value_type;
    00036         typedef const value_type &const_reference;
    00037         typedef value_type *pointer;
    00038         typedef const value_type *const_pointer;
    00039
    00040     private:
    00041         // Proxied element operations
    00042         void get_d () const {
    00043             const_pointer p = (*this) ().find_element (i_, j_);
    00044             if (p)
    00045                 d_ = *p;
    00046             else
    00047                 d_ = value_type/*zero*/();
    00048         }
    00049
    00050         void set (const value_type &s) const {
    00051             pointer p = (*this) ().find_element (i_, j_);
    00052             if (!p)
    00053                 (*this) ().insert_element (i_, j_, s);
    00054             else
    00055                 *p = s;
    00056         }
    00057
    00058     public:
    00059         // Construction and destruction
    00060         BOOST_UBLAS_INLINE
    00061         sparse_matrix_element (matrix_type &m, size_type i, size_type j):
    00062             container_reference<matrix_type> (m), i_ (i), j_ (j) {
    00063         }
    00064         BOOST_UBLAS_INLINE
    00065         sparse_matrix_element (const sparse_matrix_element &p):
    00066             container_reference<matrix_type> (p), i_ (p.i_), j_ (p.j_) {}
    00067         BOOST_UBLAS_INLINE
    00068         ~sparse_matrix_element () {
    00069         }
    00070
    00071         // Assignment
    00072         BOOST_UBLAS_INLINE
    00073         sparse_matrix_element &operator = (const sparse_matrix_element &p) {
    00074             // Overide the implict copy assignment
    00075             p.get_d ();
    00076             set (p.d_);
    00077             return *this;
    00078         }
    00079         template<class D>
    00080         BOOST_UBLAS_INLINE
    00081         sparse_matrix_element &operator = (const D &d) {
    00082             set (d);
    00083             return *this;
    00084         }
    00085         template<class D>
    00086         BOOST_UBLAS_INLINE
    00087         sparse_matrix_element &operator += (const D &d) {
    00088             get_d ();
    00089             d_ += d;
    00090             set (d_);
    00091             return *this;
    00092         }
    00093         template<class D>
    00094         BOOST_UBLAS_INLINE
    00095         sparse_matrix_element &operator -= (const D &d) {
    00096             get_d ();
    00097             d_ -= d;
    00098             set (d_);
    00099             return *this;
    00100         }
    00101         template<class D>
    00102         BOOST_UBLAS_INLINE
    00103         sparse_matrix_element &operator *= (const D &d) {
    00104             get_d ();
    00105             d_ *= d;
    00106             set (d_);
    00107             return *this;
    00108         }
    00109         template<class D>
    00110         BOOST_UBLAS_INLINE
    00111         sparse_matrix_element &operator /= (const D &d) {
    00112             get_d ();
    00113             d_ /= d;
    00114             set (d_);
    00115             return *this;
    00116         }
    00117
    00118         // Comparison
    00119         template<class D>
    00120         BOOST_UBLAS_INLINE
    00121         bool operator == (const D &d) const {
    00122             get_d ();
    00123             return d_ == d;
    00124         }
    00125         template<class D>
    00126         BOOST_UBLAS_INLINE
    00127         bool operator != (const D &d) const {
    00128             get_d ();
    00129             return d_ != d;
    00130         }
    00131
    00132         // Conversion - weak link in proxy as d_ is not a perfect alias for the element
    00133         BOOST_UBLAS_INLINE
    00134         operator const_reference () const {
    00135             get_d ();
    00136             return d_;
    00137         }
    00138
    00139         // Conversion to reference - may be invalidated
    00140         BOOST_UBLAS_INLINE
    00141         value_type& ref () const {
    00142             const pointer p = (*this) ().find_element (i_, j_);
    00143             if (!p)
    00144                 return (*this) ().insert_element (i_, j_, value_type/*zero*/());
    00145             else
    00146                 return *p;
    00147         }
    00148
    00149     private:
    00150         size_type i_;
    00151         size_type j_;
    00152         mutable value_type d_;
    00153     };
    00154
    00155     /*
    00156      * Generalise explicit reference access
    00157      */
    00158     namespace detail {
    00159         template <class V>
    00160         struct element_reference<sparse_matrix_element<V> > {
    00161             typedef typename V::value_type& reference;
    00162             static reference get_reference (const sparse_matrix_element<V>& sve)
    00163             {
    00164                 return sve.ref ();
    00165             }
    00166         };
    00167     }
    00168
    00169
    00170     template<class M>
    00171     struct type_traits<sparse_matrix_element<M> > {
    00172         typedef typename M::value_type element_type;
    00173         typedef type_traits<sparse_matrix_element<M> > self_type;
    00174         typedef typename type_traits<element_type>::value_type value_type;
    00175         typedef typename type_traits<element_type>::const_reference const_reference;
    00176         typedef sparse_matrix_element<M> reference;
    00177         typedef typename type_traits<element_type>::real_type real_type;
    00178         typedef typename type_traits<element_type>::precision_type precision_type;
    00179
    00180         static const unsigned plus_complexity = type_traits<element_type>::plus_complexity;
    00181         static const unsigned multiplies_complexity = type_traits<element_type>::multiplies_complexity;
    00182
    00183         static
    00184         BOOST_UBLAS_INLINE
    00185         real_type real (const_reference t) {
    00186             return type_traits<element_type>::real (t);
    00187         }
    00188         static
    00189         BOOST_UBLAS_INLINE
    00190         real_type imag (const_reference t) {
    00191             return type_traits<element_type>::imag (t);
    00192         }
    00193         static
    00194         BOOST_UBLAS_INLINE
    00195         value_type conj (const_reference t) {
    00196             return type_traits<element_type>::conj (t);
    00197         }
    00198
    00199         static
    00200         BOOST_UBLAS_INLINE
    00201         real_type type_abs (const_reference t) {
    00202             return type_traits<element_type>::type_abs (t);
    00203         }
    00204         static
    00205         BOOST_UBLAS_INLINE
    00206         value_type type_sqrt (const_reference t) {
    00207             return type_traits<element_type>::type_sqrt (t);
    00208         }
    00209
    00210         static
    00211         BOOST_UBLAS_INLINE
    00212         real_type norm_1 (const_reference t) {
    00213             return type_traits<element_type>::norm_1 (t);
    00214         }
    00215         static
    00216         BOOST_UBLAS_INLINE
    00217         real_type norm_2 (const_reference t) {
    00218             return type_traits<element_type>::norm_2 (t);
    00219         }
    00220         static
    00221         BOOST_UBLAS_INLINE
    00222         real_type norm_inf (const_reference t) {
    00223             return type_traits<element_type>::norm_inf (t);
    00224         }
    00225
    00226         static
    00227         BOOST_UBLAS_INLINE
    00228         bool equals (const_reference t1, const_reference t2) {
    00229             return type_traits<element_type>::equals (t1, t2);
    00230         }
    00231     };
    00232
    00233     template<class M1, class T2>
    00234     struct promote_traits<sparse_matrix_element<M1>, T2> {
    00235         typedef typename promote_traits<typename sparse_matrix_element<M1>::value_type, T2>::promote_type promote_type;
    00236     };
    00237     template<class T1, class M2>
    00238     struct promote_traits<T1, sparse_matrix_element<M2> > {
    00239         typedef typename promote_traits<T1, typename sparse_matrix_element<M2>::value_type>::promote_type promote_type;
    00240     };
    00241     template<class M1, class M2>
    00242     struct promote_traits<sparse_matrix_element<M1>, sparse_matrix_element<M2> > {
    00243         typedef typename promote_traits<typename sparse_matrix_element<M1>::value_type,
    00244                                         typename sparse_matrix_element<M2>::value_type>::promote_type promote_type;
    00245     };
    00246
    00247 #endif
    00248
    00267     template<class T, class L, class A>
    00268     class mapped_matrix:
    00269         public matrix_container<mapped_matrix<T, L, A> > {
    00270
    00271         typedef T &true_reference;
    00272         typedef T *pointer;
    00273         typedef const T * const_pointer;
    00274         typedef L layout_type;
    00275         typedef mapped_matrix<T, L, A> self_type;
    00276     public:
    00277 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
    00278         using matrix_container<self_type>::operator ();
    00279 #endif
    00280         typedef typename A::size_type size_type;
    00281         typedef typename A::difference_type difference_type;
    00282         typedef T value_type;
    00283         typedef A array_type;
    00284         typedef const T &const_reference;
    00285 #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
    00286         typedef typename detail::map_traits<A, T>::reference reference;
    00287 #else
    00288         typedef sparse_matrix_element<self_type> reference;
    00289 #endif
    00290         typedef const matrix_reference<const self_type> const_closure_type;
    00291         typedef matrix_reference<self_type> closure_type;
    00292         typedef mapped_vector<T, A> vector_temporary_type;
    00293         typedef self_type matrix_temporary_type;
    00294         typedef sparse_tag storage_category;
    00295         typedef typename L::orientation_category orientation_category;
    00296
    00297         // Construction and destruction
    00298         BOOST_UBLAS_INLINE
    00299         mapped_matrix ():
    00300             matrix_container<self_type> (),
    00301             size1_ (0), size2_ (0), data_ () {}
    00302         BOOST_UBLAS_INLINE
    00303         mapped_matrix (size_type size1, size_type size2, size_type non_zeros = 0):
    00304             matrix_container<self_type> (),
    00305             size1_ (size1), size2_ (size2), data_ () {
    00306             detail::map_reserve (data (), restrict_capacity (non_zeros));
    00307         }
    00308         BOOST_UBLAS_INLINE
    00309         mapped_matrix (const mapped_matrix &m):
    00310             matrix_container<self_type> (),
    00311             size1_ (m.size1_), size2_ (m.size2_), data_ (m.data_) {}
    00312         template<class AE>
    00313         BOOST_UBLAS_INLINE
    00314         mapped_matrix (const matrix_expression<AE> &ae, size_type non_zeros = 0):
    00315             matrix_container<self_type> (),
    00316             size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), data_ () {
    00317             detail::map_reserve (data (), restrict_capacity (non_zeros));
    00318             matrix_assign<scalar_assign> (*this, ae);
    00319         }
    00320
    00321         // Accessors
    00322         BOOST_UBLAS_INLINE
    00323         size_type size1 () const {
    00324             return size1_;
    00325         }
    00326         BOOST_UBLAS_INLINE
    00327         size_type size2 () const {
    00328             return size2_;
    00329         }
    00330         BOOST_UBLAS_INLINE
    00331         size_type nnz_capacity () const {
    00332             return detail::map_capacity (data ());
    00333         }
    00334         BOOST_UBLAS_INLINE
    00335         size_type nnz () const {
    00336             return data (). size ();
    00337         }
    00338
    00339         // Storage accessors
    00340         BOOST_UBLAS_INLINE
    00341         const array_type &data () const {
    00342             return data_;
    00343         }
    00344         BOOST_UBLAS_INLINE
    00345         array_type &data () {
    00346             return data_;
    00347         }
    00348
    00349         // Resizing
    00350     private:
    00351         BOOST_UBLAS_INLINE
    00352         size_type restrict_capacity (size_type non_zeros) const {
    00353             // Guarding against overflow - thanks to Alexei Novakov for the hint.
    00354             // non_zeros = (std::min) (non_zeros, size1_ * size2_);
    00355             if (size1_ > 0 && non_zeros / size1_ >= size2_)
    00356                 non_zeros = size1_ * size2_;
    00357             return non_zeros;
    00358         }
    00359     public:
    00360         BOOST_UBLAS_INLINE
    00361         void resize (size_type size1, size_type size2, bool preserve = true) {
    00362             // FIXME preserve unimplemented
    00363             BOOST_UBLAS_CHECK (!preserve, internal_logic ());
    00364             size1_ = size1;
    00365             size2_ = size2;
    00366             data ().clear ();
    00367         }
    00368
    00369         // Reserving
    00370         BOOST_UBLAS_INLINE
    00371         void reserve (size_type non_zeros, bool preserve = true) {
    00372             detail::map_reserve (data (), restrict_capacity (non_zeros));
    00373         }
    00374
    00375         // Element support
    00376         BOOST_UBLAS_INLINE
    00377         pointer find_element (size_type i, size_type j) {
    00378             return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i, j));
    00379         }
    00380         BOOST_UBLAS_INLINE
    00381         const_pointer find_element (size_type i, size_type j) const {
    00382             const size_type element = layout_type::element (i, size1_, j, size2_);
    00383             const_subiterator_type it (data ().find (element));
    00384             if (it == data ().end ())
    00385                 return 0;
    00386             BOOST_UBLAS_CHECK ((*it).first == element, internal_logic ());   // broken map
    00387             return &(*it).second;
    00388         }
    00389
    00390         // Element access
    00391         BOOST_UBLAS_INLINE
    00392         const_reference operator () (size_type i, size_type j) const {
    00393             const size_type element = layout_type::element (i, size1_, j, size2_);
    00394             const_subiterator_type it (data ().find (element));
    00395             if (it == data ().end ())
    00396                 return zero_;
    00397             BOOST_UBLAS_CHECK ((*it).first == element, internal_logic ());   // broken map
    00398             return (*it).second;
    00399         }
    00400         BOOST_UBLAS_INLINE
    00401         reference operator () (size_type i, size_type j) {
    00402 #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
    00403             const size_type element = layout_type::element (i, size1_, j, size2_);
    00404             std::pair<subiterator_type, bool> ii (data ().insert (typename array_type::value_type (element, value_type/*zero*/())));
    00405             BOOST_UBLAS_CHECK ((ii.first)->first == element, internal_logic ());   // broken map
    00406             return (ii.first)->second;
    00407 #else
    00408             return reference (*this, i, j);
    00409 #endif
    00410         }
    00411
    00412         // Element assingment
    00413         BOOST_UBLAS_INLINE
    00414         true_reference insert_element (size_type i, size_type j, const_reference t) {
    00415             BOOST_UBLAS_CHECK (!find_element (i, j), bad_index ());        // duplicate element
    00416             const size_type element = layout_type::element (i, size1_, j, size2_);
    00417             std::pair<subiterator_type, bool> ii (data ().insert (typename array_type::value_type (element, t)));
    00418             BOOST_UBLAS_CHECK ((ii.first)->first == element, internal_logic ());   // broken map
    00419             if (!ii.second)     // existing element
    00420                 (ii.first)->second = t;
    00421             return (ii.first)->second;
    00422         }
    00423         BOOST_UBLAS_INLINE
    00424         void erase_element (size_type i, size_type j) {
    00425             subiterator_type it = data ().find (layout_type::element (i, size1_, j, size2_));
    00426             if (it == data ().end ())
    00427                 return;
    00428             data ().erase (it);
    00429         }
    00430
    00431         // Zeroing
    00432         BOOST_UBLAS_INLINE
    00433         void clear () {
    00434             data ().clear ();
    00435         }
    00436
    00437         // Assignment
    00438         BOOST_UBLAS_INLINE
    00439         mapped_matrix &operator = (const mapped_matrix &m) {
    00440             if (this != &m) {
    00441                 size1_ = m.size1_;
    00442                 size2_ = m.size2_;
    00443                 data () = m.data ();
    00444             }
    00445             return *this;
    00446         }
    00447         template<class C>          // Container assignment without temporary
    00448         BOOST_UBLAS_INLINE
    00449         mapped_matrix &operator = (const matrix_container<C> &m) {
    00450             resize (m ().size1 (), m ().size2 (), false);
    00451             assign (m);
    00452             return *this;
    00453         }
    00454         BOOST_UBLAS_INLINE
    00455         mapped_matrix &assign_temporary (mapped_matrix &m) {
    00456             swap (m);
    00457             return *this;
    00458         }
    00459         template<class AE>
    00460         BOOST_UBLAS_INLINE
    00461         mapped_matrix &operator = (const matrix_expression<AE> &ae) {
    00462             self_type temporary (ae, detail::map_capacity (data ()));
    00463             return assign_temporary (temporary);
    00464         }
    00465         template<class AE>
    00466         BOOST_UBLAS_INLINE
    00467         mapped_matrix &assign (const matrix_expression<AE> &ae) {
    00468             matrix_assign<scalar_assign> (*this, ae);
    00469             return *this;
    00470         }
    00471         template<class AE>
    00472         BOOST_UBLAS_INLINE
    00473         mapped_matrix& operator += (const matrix_expression<AE> &ae) {
    00474             self_type temporary (*this + ae, detail::map_capacity (data ()));
    00475             return assign_temporary (temporary);
    00476         }
    00477         template<class C>          // Container assignment without temporary
    00478         BOOST_UBLAS_INLINE
    00479         mapped_matrix &operator += (const matrix_container<C> &m) {
    00480             plus_assign (m);
    00481             return *this;
    00482         }
    00483         template<class AE>
    00484         BOOST_UBLAS_INLINE
    00485         mapped_matrix &plus_assign (const matrix_expression<AE> &ae) {
    00486             matrix_assign<scalar_plus_assign> (*this, ae);
    00487             return *this;
    00488         }
    00489         template<class AE>
    00490         BOOST_UBLAS_INLINE
    00491         mapped_matrix& operator -= (const matrix_expression<AE> &ae) {
    00492             self_type temporary (*this - ae, detail::map_capacity (data ()));
    00493             return assign_temporary (temporary);
    00494         }
    00495         template<class C>          // Container assignment without temporary
    00496         BOOST_UBLAS_INLINE
    00497         mapped_matrix &operator -= (const matrix_container<C> &m) {
    00498             minus_assign (m);
    00499             return *this;
    00500         }
    00501         template<class AE>
    00502         BOOST_UBLAS_INLINE
    00503         mapped_matrix &minus_assign (const matrix_expression<AE> &ae) {
    00504             matrix_assign<scalar_minus_assign> (*this, ae);
    00505             return *this;
    00506         }
    00507         template<class AT>
    00508         BOOST_UBLAS_INLINE
    00509         mapped_matrix& operator *= (const AT &at) {
    00510             matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
    00511             return *this;
    00512         }
    00513         template<class AT>
    00514         BOOST_UBLAS_INLINE
    00515         mapped_matrix& operator /= (const AT &at) {
    00516             matrix_assign_scalar<scalar_divides_assign> (*this, at);
    00517             return *this;
    00518         }
    00519
    00520         // Swapping
    00521         BOOST_UBLAS_INLINE
    00522         void swap (mapped_matrix &m) {
    00523             if (this != &m) {
    00524                 std::swap (size1_, m.size1_);
    00525                 std::swap (size2_, m.size2_);
    00526                 data ().swap (m.data ());
    00527             }
    00528         }
    00529         BOOST_UBLAS_INLINE
    00530         friend void swap (mapped_matrix &m1, mapped_matrix &m2) {
    00531             m1.swap (m2);
    00532         }
    00533
    00534         // Iterator types
    00535     private:
    00536         // Use storage iterator
    00537         typedef typename A::const_iterator const_subiterator_type;
    00538         typedef typename A::iterator subiterator_type;
    00539
    00540         BOOST_UBLAS_INLINE
    00541         true_reference at_element (size_type i, size_type j) {
    00542             const size_type element = layout_type::element (i, size1_, j, size2_);
    00543             subiterator_type it (data ().find (element));
    00544             BOOST_UBLAS_CHECK (it != data ().end(), bad_index ());
    00545             BOOST_UBLAS_CHECK ((*it).first == element, internal_logic ());   // broken map
    00546             return it->second;
    00547         }
    00548
    00549     public:
    00550         class const_iterator1;
    00551         class iterator1;
    00552         class const_iterator2;
    00553         class iterator2;
    00554         typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
    00555         typedef reverse_iterator_base1<iterator1> reverse_iterator1;
    00556         typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
    00557         typedef reverse_iterator_base2<iterator2> reverse_iterator2;
    00558
    00559         // Element lookup
    00560         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    00561         const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const {
    00562             const_subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_)));
    00563             const_subiterator_type it_end (data ().end ());
    00564             size_type index1 = size_type (-1);
    00565             size_type index2 = size_type (-1);
    00566             while (rank == 1 && it != it_end) {
    00567                 index1 = layout_type::index_i ((*it).first, size1_, size2_);
    00568                 index2 = layout_type::index_j ((*it).first, size1_, size2_);
    00569                 if (direction > 0) {
    00570                     if ((index1 >= i && index2 == j) || (i >= size1_))
    00571                         break;
    00572                     ++ i;
    00573                 } else /* if (direction < 0) */ {
    00574                     if ((index1 <= i && index2 == j) || (i == 0))
    00575                         break;
    00576                     -- i;
    00577                 }
    00578                 it = data ().lower_bound (layout_type::address (i, size1_, j, size2_));
    00579             }
    00580             if (rank == 1 && index2 != j) {
    00581                 if (direction > 0)
    00582                     i = size1_;
    00583                 else /* if (direction < 0) */
    00584                     i = 0;
    00585                 rank = 0;
    00586             }
    00587             return const_iterator1 (*this, rank, i, j, it);
    00588         }
    00589         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    00590         iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) {
    00591             subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_)));
    00592             subiterator_type it_end (data ().end ());
    00593             size_type index1 = size_type (-1);
    00594             size_type index2 = size_type (-1);
    00595             while (rank == 1 && it != it_end) {
    00596                 index1 = layout_type::index_i ((*it).first, size1_, size2_);
    00597                 index2 = layout_type::index_j ((*it).first, size1_, size2_);
    00598                 if (direction > 0) {
    00599                     if ((index1 >= i && index2 == j) || (i >= size1_))
    00600                         break;
    00601                     ++ i;
    00602                 } else /* if (direction < 0) */ {
    00603                     if ((index1 <= i && index2 == j) || (i == 0))
    00604                         break;
    00605                     -- i;
    00606                 }
    00607                 it = data ().lower_bound (layout_type::address (i, size1_, j, size2_));
    00608             }
    00609             if (rank == 1 && index2 != j) {
    00610                 if (direction > 0)
    00611                     i = size1_;
    00612                 else /* if (direction < 0) */
    00613                     i = 0;
    00614                 rank = 0;
    00615             }
    00616             return iterator1 (*this, rank, i, j, it);
    00617         }
    00618         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    00619         const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const {
    00620             const_subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_)));
    00621             const_subiterator_type it_end (data ().end ());
    00622             size_type index1 = size_type (-1);
    00623             size_type index2 = size_type (-1);
    00624             while (rank == 1 && it != it_end) {
    00625                 index1 = layout_type::index_i ((*it).first, size1_, size2_);
    00626                 index2 = layout_type::index_j ((*it).first, size1_, size2_);
    00627                 if (direction > 0) {
    00628                     if ((index2 >= j && index1 == i) || (j >= size2_))
    00629                         break;
    00630                     ++ j;
    00631                 } else /* if (direction < 0) */ {
    00632                     if ((index2 <= j && index1 == i) || (j == 0))
    00633                         break;
    00634                     -- j;
    00635                 }
    00636                 it = data ().lower_bound (layout_type::address (i, size1_, j, size2_));
    00637             }
    00638             if (rank == 1 && index1 != i) {
    00639                 if (direction > 0)
    00640                     j = size2_;
    00641                 else /* if (direction < 0) */
    00642                     j = 0;
    00643                 rank = 0;
    00644             }
    00645             return const_iterator2 (*this, rank, i, j, it);
    00646         }
    00647         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    00648         iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) {
    00649             subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_)));
    00650             subiterator_type it_end (data ().end ());
    00651             size_type index1 = size_type (-1);
    00652             size_type index2 = size_type (-1);
    00653             while (rank == 1 && it != it_end) {
    00654                 index1 = layout_type::index_i ((*it).first, size1_, size2_);
    00655                 index2 = layout_type::index_j ((*it).first, size1_, size2_);
    00656                 if (direction > 0) {
    00657                     if ((index2 >= j && index1 == i) || (j >= size2_))
    00658                         break;
    00659                     ++ j;
    00660                 } else /* if (direction < 0) */ {
    00661                     if ((index2 <= j && index1 == i) || (j == 0))
    00662                         break;
    00663                     -- j;
    00664                 }
    00665                 it = data ().lower_bound (layout_type::address (i, size1_, j, size2_));
    00666             }
    00667             if (rank == 1 && index1 != i) {
    00668                 if (direction > 0)
    00669                     j = size2_;
    00670                 else /* if (direction < 0) */
    00671                     j = 0;
    00672                 rank = 0;
    00673             }
    00674             return iterator2 (*this, rank, i, j, it);
    00675         }
    00676
    00677
    00678         class const_iterator1:
    00679             public container_const_reference<mapped_matrix>,
    00680             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
    00681                                                const_iterator1, value_type> {
    00682         public:
    00683             typedef typename mapped_matrix::value_type value_type;
    00684             typedef typename mapped_matrix::difference_type difference_type;
    00685             typedef typename mapped_matrix::const_reference reference;
    00686             typedef const typename mapped_matrix::pointer pointer;
    00687
    00688             typedef const_iterator2 dual_iterator_type;
    00689             typedef const_reverse_iterator2 dual_reverse_iterator_type;
    00690
    00691             // Construction and destruction
    00692             BOOST_UBLAS_INLINE
    00693             const_iterator1 ():
    00694                 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), it_ () {}
    00695             BOOST_UBLAS_INLINE
    00696             const_iterator1 (const self_type &m, int rank, size_type i, size_type j, const const_subiterator_type &it):
    00697                 container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), it_ (it) {}
    00698             BOOST_UBLAS_INLINE
    00699             const_iterator1 (const iterator1 &it):
    00700                 container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), it_ (it.it_) {}
    00701
    00702             // Arithmetic
    00703             BOOST_UBLAS_INLINE
    00704             const_iterator1 &operator ++ () {
    00705                 if (rank_ == 1 && layout_type::fast_i ())
    00706                     ++ it_;
    00707                 else
    00708                     *this = (*this) ().find1 (rank_, index1 () + 1, j_, 1);
    00709                 return *this;
    00710             }
    00711             BOOST_UBLAS_INLINE
    00712             const_iterator1 &operator -- () {
    00713                 if (rank_ == 1 && layout_type::fast_i ())
    00714                     -- it_;
    00715                 else
    00716                     *this = (*this) ().find1 (rank_, index1 () - 1, j_, -1);
    00717                 return *this;
    00718             }
    00719
    00720             // Dereference
    00721             BOOST_UBLAS_INLINE
    00722             const_reference operator * () const {
    00723                 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
    00724                 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
    00725                 if (rank_ == 1) {
    00726                     return (*it_).second;
    00727                 } else {
    00728                     return (*this) () (i_, j_);
    00729                 }
    00730             }
    00731
    00732 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    00733             BOOST_UBLAS_INLINE
    00734 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    00735             typename self_type::
    00736 #endif
    00737             const_iterator2 begin () const {
    00738                 const self_type &m = (*this) ();
    00739                 return m.find2 (1, index1 (), 0);
    00740             }
    00741             BOOST_UBLAS_INLINE
    00742 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    00743             typename self_type::
    00744 #endif
    00745             const_iterator2 end () const {
    00746                 const self_type &m = (*this) ();
    00747                 return m.find2 (1, index1 (), m.size2 ());
    00748             }
    00749             BOOST_UBLAS_INLINE
    00750 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    00751             typename self_type::
    00752 #endif
    00753             const_reverse_iterator2 rbegin () const {
    00754                 return const_reverse_iterator2 (end ());
    00755             }
    00756             BOOST_UBLAS_INLINE
    00757 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    00758             typename self_type::
    00759 #endif
    00760             const_reverse_iterator2 rend () const {
    00761                 return const_reverse_iterator2 (begin ());
    00762             }
    00763 #endif
    00764
    00765             // Indices
    00766             BOOST_UBLAS_INLINE
    00767             size_type index1 () const {
    00768                 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
    00769                 if (rank_ == 1) {
    00770                     const self_type &m = (*this) ();
    00771                     BOOST_UBLAS_CHECK (layout_type::index_i ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size1 (), bad_index ());
    00772                     return layout_type::index_i ((*it_).first, m.size1 (), m.size2 ());
    00773                 } else {
    00774                     return i_;
    00775                 }
    00776             }
    00777             BOOST_UBLAS_INLINE
    00778             size_type index2 () const {
    00779                 if (rank_ == 1) {
    00780                     const self_type &m = (*this) ();
    00781                     BOOST_UBLAS_CHECK (layout_type::index_j ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size2 (), bad_index ());
    00782                     return layout_type::index_j ((*it_).first, m.size1 (), m.size2 ());
    00783                 } else {
    00784                     return j_;
    00785                 }
    00786             }
    00787
    00788             // Assignment
    00789             BOOST_UBLAS_INLINE
    00790             const_iterator1 &operator = (const const_iterator1 &it) {
    00791                 container_const_reference<self_type>::assign (&it ());
    00792                 rank_ = it.rank_;
    00793                 i_ = it.i_;
    00794                 j_ = it.j_;
    00795                 it_ = it.it_;
    00796                 return *this;
    00797             }
    00798
    00799             // Comparison
    00800             BOOST_UBLAS_INLINE
    00801             bool operator == (const const_iterator1 &it) const {
    00802                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
    00803                 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
    00804                 if (rank_ == 1 || it.rank_ == 1) {
    00805                     return it_ == it.it_;
    00806                 } else {
    00807                     return i_ == it.i_ && j_ == it.j_;
    00808                 }
    00809             }
    00810
    00811         private:
    00812             int rank_;
    00813             size_type i_;
    00814             size_type j_;
    00815             const_subiterator_type it_;
    00816         };
    00817
    00818         BOOST_UBLAS_INLINE
    00819         const_iterator1 begin1 () const {
    00820             return find1 (0, 0, 0);
    00821         }
    00822         BOOST_UBLAS_INLINE
    00823         const_iterator1 end1 () const {
    00824             return find1 (0, size1_, 0);
    00825         }
    00826
    00827         class iterator1:
    00828             public container_reference<mapped_matrix>,
    00829             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
    00830                                                iterator1, value_type> {
    00831         public:
    00832             typedef typename mapped_matrix::value_type value_type;
    00833             typedef typename mapped_matrix::difference_type difference_type;
    00834             typedef typename mapped_matrix::true_reference reference;
    00835             typedef typename mapped_matrix::pointer pointer;
    00836
    00837             typedef iterator2 dual_iterator_type;
    00838             typedef reverse_iterator2 dual_reverse_iterator_type;
    00839
    00840             // Construction and destruction
    00841             BOOST_UBLAS_INLINE
    00842             iterator1 ():
    00843                 container_reference<self_type> (), rank_ (), i_ (), j_ (), it_ () {}
    00844             BOOST_UBLAS_INLINE
    00845             iterator1 (self_type &m, int rank, size_type i, size_type j, const subiterator_type &it):
    00846                 container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), it_ (it) {}
    00847
    00848             // Arithmetic
    00849             BOOST_UBLAS_INLINE
    00850             iterator1 &operator ++ () {
    00851                 if (rank_ == 1 && layout_type::fast_i ())
    00852                     ++ it_;
    00853                 else
    00854                     *this = (*this) ().find1 (rank_, index1 () + 1, j_, 1);
    00855                 return *this;
    00856             }
    00857             BOOST_UBLAS_INLINE
    00858             iterator1 &operator -- () {
    00859                 if (rank_ == 1 && layout_type::fast_i ())
    00860                     -- it_;
    00861                 else
    00862                     *this = (*this) ().find1 (rank_, index1 () - 1, j_, -1);
    00863                 return *this;
    00864             }
    00865
    00866             // Dereference
    00867             BOOST_UBLAS_INLINE
    00868             reference operator * () const {
    00869                 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
    00870                 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
    00871                 if (rank_ == 1) {
    00872                     return (*it_).second;
    00873                 } else {
    00874                     return (*this) ().at_element (i_, j_);
    00875                 }
    00876             }
    00877
    00878 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    00879             BOOST_UBLAS_INLINE
    00880 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    00881             typename self_type::
    00882 #endif
    00883             iterator2 begin () const {
    00884                 self_type &m = (*this) ();
    00885                 return m.find2 (1, index1 (), 0);
    00886             }
    00887             BOOST_UBLAS_INLINE
    00888 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    00889             typename self_type::
    00890 #endif
    00891             iterator2 end () const {
    00892                 self_type &m = (*this) ();
    00893                 return m.find2 (1, index1 (), m.size2 ());
    00894             }
    00895             BOOST_UBLAS_INLINE
    00896 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    00897             typename self_type::
    00898 #endif
    00899             reverse_iterator2 rbegin () const {
    00900                 return reverse_iterator2 (end ());
    00901             }
    00902             BOOST_UBLAS_INLINE
    00903 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    00904             typename self_type::
    00905 #endif
    00906             reverse_iterator2 rend () const {
    00907                 return reverse_iterator2 (begin ());
    00908             }
    00909 #endif
    00910
    00911             // Indices
    00912             BOOST_UBLAS_INLINE
    00913             size_type index1 () const {
    00914                 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
    00915                 if (rank_ == 1) {
    00916                     const self_type &m = (*this) ();
    00917                     BOOST_UBLAS_CHECK (layout_type::index_i ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size1 (), bad_index ());
    00918                     return layout_type::index_i ((*it_).first, m.size1 (), m.size2 ());
    00919                 } else {
    00920                     return i_;
    00921                 }
    00922             }
    00923             BOOST_UBLAS_INLINE
    00924             size_type index2 () const {
    00925                 if (rank_ == 1) {
    00926                     const self_type &m = (*this) ();
    00927                     BOOST_UBLAS_CHECK (layout_type::index_j ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size2 (), bad_index ());
    00928                     return layout_type::index_j ((*it_).first, m.size1 (), m.size2 ());
    00929                 } else {
    00930                     return j_;
    00931                 }
    00932             }
    00933
    00934             // Assignment
    00935             BOOST_UBLAS_INLINE
    00936             iterator1 &operator = (const iterator1 &it) {
    00937                 container_reference<self_type>::assign (&it ());
    00938                 rank_ = it.rank_;
    00939                 i_ = it.i_;
    00940                 j_ = it.j_;
    00941                 it_ = it.it_;
    00942                 return *this;
    00943             }
    00944
    00945             // Comparison
    00946             BOOST_UBLAS_INLINE
    00947             bool operator == (const iterator1 &it) const {
    00948                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
    00949                 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
    00950                 if (rank_ == 1 || it.rank_ == 1) {
    00951                     return it_ == it.it_;
    00952                 } else {
    00953                     return i_ == it.i_ && j_ == it.j_;
    00954                 }
    00955             }
    00956
    00957         private:
    00958             int rank_;
    00959             size_type i_;
    00960             size_type j_;
    00961             subiterator_type it_;
    00962
    00963             friend class const_iterator1;
    00964         };
    00965
    00966         BOOST_UBLAS_INLINE
    00967         iterator1 begin1 () {
    00968             return find1 (0, 0, 0);
    00969         }
    00970         BOOST_UBLAS_INLINE
    00971         iterator1 end1 () {
    00972             return find1 (0, size1_, 0);
    00973         }
    00974
    00975         class const_iterator2:
    00976             public container_const_reference<mapped_matrix>,
    00977             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
    00978                                                const_iterator2, value_type> {
    00979         public:
    00980             typedef typename mapped_matrix::value_type value_type;
    00981             typedef typename mapped_matrix::difference_type difference_type;
    00982             typedef typename mapped_matrix::const_reference reference;
    00983             typedef const typename mapped_matrix::pointer pointer;
    00984
    00985             typedef const_iterator1 dual_iterator_type;
    00986             typedef const_reverse_iterator1 dual_reverse_iterator_type;
    00987
    00988             // Construction and destruction
    00989             BOOST_UBLAS_INLINE
    00990             const_iterator2 ():
    00991                 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), it_ () {}
    00992             BOOST_UBLAS_INLINE
    00993             const_iterator2 (const self_type &m, int rank, size_type i, size_type j, const const_subiterator_type &it):
    00994                 container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), it_ (it) {}
    00995             BOOST_UBLAS_INLINE
    00996             const_iterator2 (const iterator2 &it):
    00997                 container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), it_ (it.it_) {}
    00998
    00999             // Arithmetic
    01000             BOOST_UBLAS_INLINE
    01001             const_iterator2 &operator ++ () {
    01002                 if (rank_ == 1 && layout_type::fast_j ())
    01003                     ++ it_;
    01004                 else
    01005                     *this = (*this) ().find2 (rank_, i_, index2 () + 1, 1);
    01006                 return *this;
    01007             }
    01008             BOOST_UBLAS_INLINE
    01009             const_iterator2 &operator -- () {
    01010                 if (rank_ == 1 && layout_type::fast_j ())
    01011                     -- it_;
    01012                 else
    01013                     *this = (*this) ().find2 (rank_, i_, index2 () - 1, -1);
    01014                 return *this;
    01015             }
    01016
    01017             // Dereference
    01018             BOOST_UBLAS_INLINE
    01019             const_reference operator * () const {
    01020                 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
    01021                 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
    01022                 if (rank_ == 1) {
    01023                     return (*it_).second;
    01024                 } else {
    01025                     return (*this) () (i_, j_);
    01026                 }
    01027             }
    01028
    01029 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    01030             BOOST_UBLAS_INLINE
    01031 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    01032             typename self_type::
    01033 #endif
    01034             const_iterator1 begin () const {
    01035                 const self_type &m = (*this) ();
    01036                 return m.find1 (1, 0, index2 ());
    01037             }
    01038             BOOST_UBLAS_INLINE
    01039 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    01040             typename self_type::
    01041 #endif
    01042             const_iterator1 end () const {
    01043                 const self_type &m = (*this) ();
    01044                 return m.find1 (1, m.size1 (), index2 ());
    01045             }
    01046             BOOST_UBLAS_INLINE
    01047 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    01048             typename self_type::
    01049 #endif
    01050             const_reverse_iterator1 rbegin () const {
    01051                 return const_reverse_iterator1 (end ());
    01052             }
    01053             BOOST_UBLAS_INLINE
    01054 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    01055             typename self_type::
    01056 #endif
    01057             const_reverse_iterator1 rend () const {
    01058                 return const_reverse_iterator1 (begin ());
    01059             }
    01060 #endif
    01061
    01062             // Indices
    01063             BOOST_UBLAS_INLINE
    01064             size_type index1 () const {
    01065                 if (rank_ == 1) {
    01066                     const self_type &m = (*this) ();
    01067                     BOOST_UBLAS_CHECK (layout_type::index_i ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size1 (), bad_index ());
    01068                     return layout_type::index_i ((*it_).first, m.size1 (), m.size2 ());
    01069                 } else {
    01070                     return i_;
    01071                 }
    01072             }
    01073             BOOST_UBLAS_INLINE
    01074             size_type index2 () const {
    01075                 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
    01076                 if (rank_ == 1) {
    01077                     const self_type &m = (*this) ();
    01078                     BOOST_UBLAS_CHECK (layout_type::index_j ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size2 (), bad_index ());
    01079                     return layout_type::index_j ((*it_).first, m.size1 (), m.size2 ());
    01080                 } else {
    01081                     return j_;
    01082                 }
    01083             }
    01084
    01085             // Assignment
    01086             BOOST_UBLAS_INLINE
    01087             const_iterator2 &operator = (const const_iterator2 &it) {
    01088                 container_const_reference<self_type>::assign (&it ());
    01089                 rank_ = it.rank_;
    01090                 i_ = it.i_;
    01091                 j_ = it.j_;
    01092                 it_ = it.it_;
    01093                 return *this;
    01094             }
    01095
    01096             // Comparison
    01097             BOOST_UBLAS_INLINE
    01098             bool operator == (const const_iterator2 &it) const {
    01099                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
    01100                 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
    01101                 if (rank_ == 1 || it.rank_ == 1) {
    01102                     return it_ == it.it_;
    01103                 } else {
    01104                     return i_ == it.i_ && j_ == it.j_;
    01105                 }
    01106             }
    01107
    01108         private:
    01109             int rank_;
    01110             size_type i_;
    01111             size_type j_;
    01112             const_subiterator_type it_;
    01113         };
    01114
    01115         BOOST_UBLAS_INLINE
    01116         const_iterator2 begin2 () const {
    01117             return find2 (0, 0, 0);
    01118         }
    01119         BOOST_UBLAS_INLINE
    01120         const_iterator2 end2 () const {
    01121             return find2 (0, 0, size2_);
    01122         }
    01123
    01124         class iterator2:
    01125             public container_reference<mapped_matrix>,
    01126             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
    01127                                                iterator2, value_type> {
    01128         public:
    01129             typedef typename mapped_matrix::value_type value_type;
    01130             typedef typename mapped_matrix::difference_type difference_type;
    01131             typedef typename mapped_matrix::true_reference reference;
    01132             typedef typename mapped_matrix::pointer pointer;
    01133
    01134             typedef iterator1 dual_iterator_type;
    01135             typedef reverse_iterator1 dual_reverse_iterator_type;
    01136
    01137             // Construction and destruction
    01138             BOOST_UBLAS_INLINE
    01139             iterator2 ():
    01140                 container_reference<self_type> (), rank_ (), i_ (), j_ (), it_ () {}
    01141             BOOST_UBLAS_INLINE
    01142             iterator2 (self_type &m, int rank, size_type i, size_type j, const subiterator_type &it):
    01143                 container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), it_ (it) {}
    01144
    01145             // Arithmetic
    01146             BOOST_UBLAS_INLINE
    01147             iterator2 &operator ++ () {
    01148                 if (rank_ == 1 && layout_type::fast_j ())
    01149                     ++ it_;
    01150                 else
    01151                     *this = (*this) ().find2 (rank_, i_, index2 () + 1, 1);
    01152                 return *this;
    01153             }
    01154             BOOST_UBLAS_INLINE
    01155             iterator2 &operator -- () {
    01156                 if (rank_ == 1 && layout_type::fast_j ())
    01157                     -- it_;
    01158                 else
    01159                     *this = (*this) ().find2 (rank_, i_, index2 () - 1, -1);
    01160                 return *this;
    01161             }
    01162
    01163             // Dereference
    01164             BOOST_UBLAS_INLINE
    01165             reference operator * () const {
    01166                 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
    01167                 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
    01168                 if (rank_ == 1) {
    01169                     return (*it_).second;
    01170                 } else {
    01171                     return (*this) ().at_element (i_, j_);
    01172                 }
    01173             }
    01174
    01175 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    01176             BOOST_UBLAS_INLINE
    01177 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    01178             typename self_type::
    01179 #endif
    01180             iterator1 begin () const {
    01181                 self_type &m = (*this) ();
    01182                 return m.find1 (1, 0, index2 ());
    01183             }
    01184             BOOST_UBLAS_INLINE
    01185 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    01186             typename self_type::
    01187 #endif
    01188             iterator1 end () const {
    01189                 self_type &m = (*this) ();
    01190                 return m.find1 (1, m.size1 (), index2 ());
    01191             }
    01192             BOOST_UBLAS_INLINE
    01193 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    01194             typename self_type::
    01195 #endif
    01196             reverse_iterator1 rbegin () const {
    01197                 return reverse_iterator1 (end ());
    01198             }
    01199             BOOST_UBLAS_INLINE
    01200 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    01201             typename self_type::
    01202 #endif
    01203             reverse_iterator1 rend () const {
    01204                 return reverse_iterator1 (begin ());
    01205             }
    01206 #endif
    01207
    01208             // Indices
    01209             BOOST_UBLAS_INLINE
    01210             size_type index1 () const {
    01211                 if (rank_ == 1) {
    01212                     const self_type &m = (*this) ();
    01213                     BOOST_UBLAS_CHECK (layout_type::index_i ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size1 (), bad_index ());
    01214                     return layout_type::index_i ((*it_).first, m.size1 (), m.size2 ());
    01215                 } else {
    01216                     return i_;
    01217                 }
    01218             }
    01219             BOOST_UBLAS_INLINE
    01220             size_type index2 () const {
    01221                 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
    01222                 if (rank_ == 1) {
    01223                     const self_type &m = (*this) ();
    01224                     BOOST_UBLAS_CHECK (layout_type::index_j ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size2 (), bad_index ());
    01225                     return layout_type::index_j ((*it_).first, m.size1 (), m.size2 ());
    01226                 } else {
    01227                     return j_;
    01228                 }
    01229             }
    01230
    01231             // Assignment
    01232             BOOST_UBLAS_INLINE
    01233             iterator2 &operator = (const iterator2 &it) {
    01234                 container_reference<self_type>::assign (&it ());
    01235                 rank_ = it.rank_;
    01236                 i_ = it.i_;
    01237                 j_ = it.j_;
    01238                 it_ = it.it_;
    01239                 return *this;
    01240             }
    01241
    01242             // Comparison
    01243             BOOST_UBLAS_INLINE
    01244             bool operator == (const iterator2 &it) const {
    01245                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
    01246                 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
    01247                 if (rank_ == 1 || it.rank_ == 1) {
    01248                     return it_ == it.it_;
    01249                 } else {
    01250                     return i_ == it.i_ && j_ == it.j_;
    01251                 }
    01252             }
    01253
    01254         private:
    01255             int rank_;
    01256             size_type i_;
    01257             size_type j_;
    01258             subiterator_type it_;
    01259
    01260             friend class const_iterator2;
    01261         };
    01262
    01263         BOOST_UBLAS_INLINE
    01264         iterator2 begin2 () {
    01265             return find2 (0, 0, 0);
    01266         }
    01267         BOOST_UBLAS_INLINE
    01268         iterator2 end2 () {
    01269             return find2 (0, 0, size2_);
    01270         }
    01271
    01272         // Reverse iterators
    01273
    01274         BOOST_UBLAS_INLINE
    01275         const_reverse_iterator1 rbegin1 () const {
    01276             return const_reverse_iterator1 (end1 ());
    01277         }
    01278         BOOST_UBLAS_INLINE
    01279         const_reverse_iterator1 rend1 () const {
    01280             return const_reverse_iterator1 (begin1 ());
    01281         }
    01282
    01283         BOOST_UBLAS_INLINE
    01284         reverse_iterator1 rbegin1 () {
    01285             return reverse_iterator1 (end1 ());
    01286         }
    01287         BOOST_UBLAS_INLINE
    01288         reverse_iterator1 rend1 () {
    01289             return reverse_iterator1 (begin1 ());
    01290         }
    01291
    01292         BOOST_UBLAS_INLINE
    01293         const_reverse_iterator2 rbegin2 () const {
    01294             return const_reverse_iterator2 (end2 ());
    01295         }
    01296         BOOST_UBLAS_INLINE
    01297         const_reverse_iterator2 rend2 () const {
    01298             return const_reverse_iterator2 (begin2 ());
    01299         }
    01300
    01301         BOOST_UBLAS_INLINE
    01302         reverse_iterator2 rbegin2 () {
    01303             return reverse_iterator2 (end2 ());
    01304         }
    01305         BOOST_UBLAS_INLINE
    01306         reverse_iterator2 rend2 () {
    01307             return reverse_iterator2 (begin2 ());
    01308         }
    01309
    01310          // Serialization
    01311         template<class Archive>
    01312         void serialize(Archive & ar, const unsigned int /* file_version */){
    01313             serialization::collection_size_type s1 (size1_);
    01314             serialization::collection_size_type s2 (size2_);
    01315             ar & serialization::make_nvp("size1",s1);
    01316             ar & serialization::make_nvp("size2",s2);
    01317             if (Archive::is_loading::value) {
    01318                 size1_ = s1;
    01319                 size2_ = s2;
    01320             }
    01321             ar & serialization::make_nvp("data", data_);
    01322         }
    01323
    01324     private:
    01325         size_type size1_;
    01326         size_type size2_;
    01327         array_type data_;
    01328         static const value_type zero_;
    01329     };
    01330
    01331     template<class T, class L, class A>
    01332     const typename mapped_matrix<T, L, A>::value_type mapped_matrix<T, L, A>::zero_ = value_type/*zero*/();
    01333
    01334
    01335     // Vector index map based sparse matrix class
    01336     template<class T, class L, class A>
    01337     class mapped_vector_of_mapped_vector:
    01338         public matrix_container<mapped_vector_of_mapped_vector<T, L, A> > {
    01339
    01340         typedef T &true_reference;
    01341         typedef T *pointer;
    01342         typedef const T *const_pointer;
    01343         typedef A array_type;
    01344         typedef const A const_array_type;
    01345         typedef L layout_type;
    01346         typedef mapped_vector_of_mapped_vector<T, L, A> self_type;
    01347     public:
    01348 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
    01349         using matrix_container<self_type>::operator ();
    01350 #endif
    01351         typedef typename A::size_type size_type;
    01352         typedef typename A::difference_type difference_type;
    01353         typedef T value_type;
    01354         typedef const T &const_reference;
    01355 #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
    01356         typedef typename detail::map_traits<typename A::data_value_type, T>::reference reference;
    01357 #else
    01358         typedef sparse_matrix_element<self_type> reference;
    01359 #endif
    01360         typedef const matrix_reference<const self_type> const_closure_type;
    01361         typedef matrix_reference<self_type> closure_type;
    01362         typedef mapped_vector<T, typename A::value_type> vector_temporary_type;
    01363         typedef self_type matrix_temporary_type;
    01364         typedef typename A::value_type::second_type vector_data_value_type;
    01365         typedef sparse_tag storage_category;
    01366         typedef typename L::orientation_category orientation_category;
    01367
    01368         // Construction and destruction
    01369         BOOST_UBLAS_INLINE
    01370         mapped_vector_of_mapped_vector ():
    01371             matrix_container<self_type> (),
    01372             size1_ (0), size2_ (0), data_ () {
    01373             data_ [layout_type::size_M (size1_, size2_)] = vector_data_value_type ();
    01374         }
    01375         BOOST_UBLAS_INLINE
    01376         mapped_vector_of_mapped_vector (size_type size1, size_type size2, size_type non_zeros = 0):
    01377             matrix_container<self_type> (),
    01378             size1_ (size1), size2_ (size2), data_ () {
    01379             data_ [layout_type::size_M (size1_, size2_)] = vector_data_value_type ();
    01380         }
    01381         BOOST_UBLAS_INLINE
    01382         mapped_vector_of_mapped_vector (const mapped_vector_of_mapped_vector &m):
    01383             matrix_container<self_type> (),
    01384             size1_ (m.size1_), size2_ (m.size2_), data_ (m.data_) {}
    01385         template<class AE>
    01386         BOOST_UBLAS_INLINE
    01387         mapped_vector_of_mapped_vector (const matrix_expression<AE> &ae, size_type non_zeros = 0):
    01388             matrix_container<self_type> (),
    01389             size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), data_ () {
    01390             data_ [layout_type::size_M (size1_, size2_)] = vector_data_value_type ();
    01391             matrix_assign<scalar_assign> (*this, ae);
    01392         }
    01393
    01394         // Accessors
    01395         BOOST_UBLAS_INLINE
    01396         size_type size1 () const {
    01397             return size1_;
    01398         }
    01399         BOOST_UBLAS_INLINE
    01400         size_type size2 () const {
    01401             return size2_;
    01402         }
    01403         BOOST_UBLAS_INLINE
    01404         size_type nnz_capacity () const {
    01405             size_type non_zeros = 0;
    01406             for (vector_const_subiterator_type itv = data_ ().begin (); itv != data_ ().end (); ++ itv)
    01407                 non_zeros += detail::map_capacity (*itv);
    01408             return non_zeros;
    01409         }
    01410         BOOST_UBLAS_INLINE
    01411         size_type nnz () const {
    01412             size_type filled = 0;
    01413             for (vector_const_subiterator_type itv = data_ ().begin (); itv != data_ ().end (); ++ itv)
    01414                 filled += (*itv).size ();
    01415             return filled;
    01416         }
    01417
    01418         // Storage accessors
    01419         BOOST_UBLAS_INLINE
    01420         const_array_type &data () const {
    01421             return data_;
    01422         }
    01423         BOOST_UBLAS_INLINE
    01424         array_type &data () {
    01425             return data_;
    01426         }
    01427
    01428         // Resizing
    01429         BOOST_UBLAS_INLINE
    01430         void resize (size_type size1, size_type size2, bool preserve = true) {
    01431             // FIXME preserve unimplemented
    01432             BOOST_UBLAS_CHECK (!preserve, internal_logic ());
    01433             size1_ = size1;
    01434             size2_ = size2;
    01435             data ().clear ();
    01436             data () [layout_type::size_M (size1_, size2_)] = vector_data_value_type ();
    01437         }
    01438
    01439         // Element support
    01440         BOOST_UBLAS_INLINE
    01441         pointer find_element (size_type i, size_type j) {
    01442             return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i, j));
    01443         }
    01444         BOOST_UBLAS_INLINE
    01445         const_pointer find_element (size_type i, size_type j) const {
    01446             const size_type element1 = layout_type::index_M (i, j);
    01447             const size_type element2 = layout_type::index_m (i, j);
    01448             vector_const_subiterator_type itv (data ().find (element1));
    01449             if (itv == data ().end ())
    01450                 return 0;
    01451             BOOST_UBLAS_CHECK ((*itv).first == element1, internal_logic ());   // broken map
    01452             const_subiterator_type it ((*itv).second.find (element2));
    01453             if (it == (*itv).second.end ())
    01454                 return 0;
    01455             BOOST_UBLAS_CHECK ((*it).first == element2, internal_logic ());   // broken map
    01456             return &(*it).second;
    01457         }
    01458
    01459         // Element access
    01460         BOOST_UBLAS_INLINE
    01461         const_reference operator () (size_type i, size_type j) const {
    01462             const size_type element1 = layout_type::index_M (i, j);
    01463             const size_type element2 = layout_type::index_m (i, j);
    01464             vector_const_subiterator_type itv (data ().find (element1));
    01465             if (itv == data ().end ())
    01466                 return zero_;
    01467             BOOST_UBLAS_CHECK ((*itv).first == element1, internal_logic ());   // broken map
    01468             const_subiterator_type it ((*itv).second.find (element2));
    01469             if (it == (*itv).second.end ())
    01470                 return zero_;
    01471             BOOST_UBLAS_CHECK ((*itv).first == element1, internal_logic ());   // broken map
    01472             return (*it).second;
    01473         }
    01474         BOOST_UBLAS_INLINE
    01475         reference operator () (size_type i, size_type j) {
    01476 #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
    01477             const size_type element1 = layout_type::index_M (i, j);
    01478             const size_type element2 = layout_type::index_m (i, j);
    01479             vector_data_value_type& vd (data () [element1]);
    01480             std::pair<subiterator_type, bool> ii (vd.insert (typename array_type::value_type::second_type::value_type (element2, value_type/*zero*/())));
    01481             BOOST_UBLAS_CHECK ((ii.first)->first == element2, internal_logic ());   // broken map
    01482             return (ii.first)->second;
    01483 #else
    01484             return reference (*this, i, j);
    01485 #endif
    01486         }
    01487
    01488         // Element assignment
    01489         BOOST_UBLAS_INLINE
    01490         true_reference insert_element (size_type i, size_type j, const_reference t) {
    01491             BOOST_UBLAS_CHECK (!find_element (i, j), bad_index ());        // duplicate element
    01492             const size_type element1 = layout_type::index_M (i, j);
    01493             const size_type element2 = layout_type::index_m (i, j);
    01494
    01495             vector_data_value_type& vd (data () [element1]);
    01496             std::pair<subiterator_type, bool> ii (vd.insert (typename vector_data_value_type::value_type (element2, t)));
    01497             BOOST_UBLAS_CHECK ((ii.first)->first == element2, internal_logic ());   // broken map
    01498             if (!ii.second)     // existing element
    01499                 (ii.first)->second = t;
    01500             return (ii.first)->second;
    01501         }
    01502         BOOST_UBLAS_INLINE
    01503         void erase_element (size_type i, size_type j) {
    01504             vector_subiterator_type itv (data ().find (layout_type::index_M (i, j)));
    01505             if (itv == data ().end ())
    01506                 return;
    01507             subiterator_type it ((*itv).second.find (layout_type::index_m (i, j)));
    01508             if (it == (*itv).second.end ())
    01509                 return;
    01510             (*itv).second.erase (it);
    01511         }
    01512
    01513         // Zeroing
    01514         BOOST_UBLAS_INLINE
    01515         void clear () {
    01516             data ().clear ();
    01517             data_ [layout_type::size_M (size1_, size2_)] = vector_data_value_type ();
    01518         }
    01519
    01520         // Assignment
    01521         BOOST_UBLAS_INLINE
    01522         mapped_vector_of_mapped_vector &operator = (const mapped_vector_of_mapped_vector &m) {
    01523             if (this != &m) {
    01524                 size1_ = m.size1_;
    01525                 size2_ = m.size2_;
    01526                 data () = m.data ();
    01527             }
    01528             return *this;
    01529         }
    01530         template<class C>          // Container assignment without temporary
    01531         BOOST_UBLAS_INLINE
    01532         mapped_vector_of_mapped_vector &operator = (const matrix_container<C> &m) {
    01533             resize (m ().size1 (), m ().size2 ());
    01534             assign (m);
    01535             return *this;
    01536         }
    01537         BOOST_UBLAS_INLINE
    01538         mapped_vector_of_mapped_vector &assign_temporary (mapped_vector_of_mapped_vector &m) {
    01539             swap (m);
    01540             return *this;
    01541         }
    01542         template<class AE>
    01543         BOOST_UBLAS_INLINE
    01544         mapped_vector_of_mapped_vector &operator = (const matrix_expression<AE> &ae) {
    01545             self_type temporary (ae);
    01546             return assign_temporary (temporary);
    01547         }
    01548         template<class AE>
    01549         BOOST_UBLAS_INLINE
    01550         mapped_vector_of_mapped_vector &assign (const matrix_expression<AE> &ae) {
    01551             matrix_assign<scalar_assign> (*this, ae);
    01552             return *this;
    01553         }
    01554         template<class AE>
    01555         BOOST_UBLAS_INLINE
    01556         mapped_vector_of_mapped_vector& operator += (const matrix_expression<AE> &ae) {
    01557             self_type temporary (*this + ae);
    01558             return assign_temporary (temporary);
    01559         }
    01560         template<class C>          // Container assignment without temporary
    01561         BOOST_UBLAS_INLINE
    01562         mapped_vector_of_mapped_vector &operator += (const matrix_container<C> &m) {
    01563             plus_assign (m);
    01564             return *this;
    01565         }
    01566         template<class AE>
    01567         BOOST_UBLAS_INLINE
    01568         mapped_vector_of_mapped_vector &plus_assign (const matrix_expression<AE> &ae) {
    01569             matrix_assign<scalar_plus_assign> (*this, ae);
    01570             return *this;
    01571         }
    01572         template<class AE>
    01573         BOOST_UBLAS_INLINE
    01574         mapped_vector_of_mapped_vector& operator -= (const matrix_expression<AE> &ae) {
    01575             self_type temporary (*this - ae);
    01576             return assign_temporary (temporary);
    01577         }
    01578         template<class C>          // Container assignment without temporary
    01579         BOOST_UBLAS_INLINE
    01580         mapped_vector_of_mapped_vector &operator -= (const matrix_container<C> &m) {
    01581             minus_assign (m);
    01582             return *this;
    01583         }
    01584         template<class AE>
    01585         BOOST_UBLAS_INLINE
    01586         mapped_vector_of_mapped_vector &minus_assign (const matrix_expression<AE> &ae) {
    01587             matrix_assign<scalar_minus_assign> (*this, ae);
    01588             return *this;
    01589         }
    01590         template<class AT>
    01591         BOOST_UBLAS_INLINE
    01592         mapped_vector_of_mapped_vector& operator *= (const AT &at) {
    01593             matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
    01594             return *this;
    01595         }
    01596         template<class AT>
    01597         BOOST_UBLAS_INLINE
    01598         mapped_vector_of_mapped_vector& operator /= (const AT &at) {
    01599             matrix_assign_scalar<scalar_divides_assign> (*this, at);
    01600             return *this;
    01601         }
    01602
    01603         // Swapping
    01604         BOOST_UBLAS_INLINE
    01605         void swap (mapped_vector_of_mapped_vector &m) {
    01606             if (this != &m) {
    01607                 std::swap (size1_, m.size1_);
    01608                 std::swap (size2_, m.size2_);
    01609                 data ().swap (m.data ());
    01610             }
    01611         }
    01612         BOOST_UBLAS_INLINE
    01613         friend void swap (mapped_vector_of_mapped_vector &m1, mapped_vector_of_mapped_vector &m2) {
    01614             m1.swap (m2);
    01615         }
    01616
    01617         // Iterator types
    01618     private:
    01619         // Use storage iterators
    01620         typedef typename A::const_iterator vector_const_subiterator_type;
    01621         typedef typename A::iterator vector_subiterator_type;
    01622         typedef typename A::value_type::second_type::const_iterator const_subiterator_type;
    01623         typedef typename A::value_type::second_type::iterator subiterator_type;
    01624
    01625         BOOST_UBLAS_INLINE
    01626         true_reference at_element (size_type i, size_type j) {
    01627             const size_type element1 = layout_type::index_M (i, j);
    01628             const size_type element2 = layout_type::index_m (i, j);
    01629             vector_subiterator_type itv (data ().find (element1));
    01630             BOOST_UBLAS_CHECK (itv != data ().end(), bad_index ());
    01631             BOOST_UBLAS_CHECK ((*itv).first == element1, internal_logic ());   // broken map
    01632             subiterator_type it ((*itv).second.find (element2));
    01633             BOOST_UBLAS_CHECK (it != (*itv).second.end (), bad_index ());
    01634             BOOST_UBLAS_CHECK ((*it).first == element2, internal_logic ());    // broken map
    01635
    01636             return it->second;
    01637         }
    01638
    01639     public:
    01640         class const_iterator1;
    01641         class iterator1;
    01642         class const_iterator2;
    01643         class iterator2;
    01644         typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
    01645         typedef reverse_iterator_base1<iterator1> reverse_iterator1;
    01646         typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
    01647         typedef reverse_iterator_base2<iterator2> reverse_iterator2;
    01648
    01649         // Element lookup
    01650         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    01651         const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const {
    01652             BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ());
    01653             for (;;) {
    01654                 vector_const_subiterator_type itv (data ().lower_bound (layout_type::index_M (i, j)));
    01655                 vector_const_subiterator_type itv_end (data ().end ());
    01656                 if (itv == itv_end)
    01657                     return const_iterator1 (*this, rank, i, j, itv_end, (*(-- itv)).second.end ());
    01658
    01659                 const_subiterator_type it ((*itv).second.lower_bound (layout_type::index_m (i, j)));
    01660                 const_subiterator_type it_end ((*itv).second.end ());
    01661                 if (rank == 0) {
    01662                     // advance to the first available major index
    01663                     size_type M = itv->first;
    01664                     size_type m;
    01665                     if (it != it_end) {
    01666                         m = it->first;
    01667                     } else {
    01668                         m = layout_type::size_m(size1_, size2_);
    01669                     }
    01670                     size_type first_i = layout_type::index_M(M,m);
    01671                     return const_iterator1 (*this, rank, first_i, j, itv, it);
    01672                 }
    01673                 if (it != it_end && (*it).first == layout_type::index_m (i, j))
    01674                     return const_iterator1 (*this, rank, i, j, itv, it);
    01675                 if (direction > 0) {
    01676                     if (layout_type::fast_i ()) {
    01677                         if (it == it_end)
    01678                             return const_iterator1 (*this, rank, i, j, itv, it);
    01679                         i = (*it).first;
    01680                     } else {
    01681                         if (i >= size1_)
    01682                             return const_iterator1 (*this, rank, i, j, itv, it);
    01683                         ++ i;
    01684                     }
    01685                 } else /* if (direction < 0)  */ {
    01686                     if (layout_type::fast_i ()) {
    01687                         if (it == (*itv).second.begin ())
    01688                             return const_iterator1 (*this, rank, i, j, itv, it);
    01689                         -- it;
    01690                         i = (*it).first;
    01691                     } else {
    01692                         if (i == 0)
    01693                             return const_iterator1 (*this, rank, i, j, itv, it);
    01694                         -- i;
    01695                     }
    01696                 }
    01697             }
    01698         }
    01699         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    01700         iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) {
    01701             BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ());
    01702             for (;;) {
    01703                 vector_subiterator_type itv (data ().lower_bound (layout_type::index_M (i, j)));
    01704                 vector_subiterator_type itv_end (data ().end ());
    01705                 if (itv == itv_end)
    01706                     return iterator1 (*this, rank, i, j, itv_end, (*(-- itv)).second.end ());
    01707
    01708                 subiterator_type it ((*itv).second.lower_bound (layout_type::index_m (i, j)));
    01709                 subiterator_type it_end ((*itv).second.end ());
    01710                 if (rank == 0) {
    01711                     // advance to the first available major index
    01712                     size_type M = itv->first;
    01713                     size_type m;
    01714                     if (it != it_end) {
    01715                         m = it->first;
    01716                     } else {
    01717                         m = layout_type::size_m(size1_, size2_);
    01718                     }
    01719                     size_type first_i = layout_type::index_M(M,m);
    01720                     return iterator1 (*this, rank, first_i, j, itv, it);
    01721                 }
    01722                 if (it != it_end && (*it).first == layout_type::index_m (i, j))
    01723                     return iterator1 (*this, rank, i, j, itv, it);
    01724                 if (direction > 0) {
    01725                     if (layout_type::fast_i ()) {
    01726                         if (it == it_end)
    01727                             return iterator1 (*this, rank, i, j, itv, it);
    01728                         i = (*it).first;
    01729                     } else {
    01730                         if (i >= size1_)
    01731                             return iterator1 (*this, rank, i, j, itv, it);
    01732                         ++ i;
    01733                     }
    01734                 } else /* if (direction < 0)  */ {
    01735                     if (layout_type::fast_i ()) {
    01736                         if (it == (*itv).second.begin ())
    01737                             return iterator1 (*this, rank, i, j, itv, it);
    01738                         -- it;
    01739                         i = (*it).first;
    01740                     } else {
    01741                         if (i == 0)
    01742                             return iterator1 (*this, rank, i, j, itv, it);
    01743                         -- i;
    01744                     }
    01745                 }
    01746             }
    01747         }
    01748         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    01749         const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const {
    01750             BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ());
    01751             for (;;) {
    01752                 vector_const_subiterator_type itv (data ().lower_bound (layout_type::index_M (i, j)));
    01753                 vector_const_subiterator_type itv_end (data ().end ());
    01754                 if (itv == itv_end)
    01755                     return const_iterator2 (*this, rank, i, j, itv_end, (*(-- itv)).second.end ());
    01756
    01757                 const_subiterator_type it ((*itv).second.lower_bound (layout_type::index_m (i, j)));
    01758                 const_subiterator_type it_end ((*itv).second.end ());
    01759                 if (rank == 0) {
    01760                     // advance to the first available major index
    01761                     size_type M = itv->first;
    01762                     size_type m;
    01763                     if (it != it_end) {
    01764                         m = it->first;
    01765                     } else {
    01766                         m = layout_type::size_m(size1_, size2_);
    01767                     }
    01768                     size_type first_j = layout_type::index_m(M,m);
    01769                     return const_iterator2 (*this, rank, i, first_j, itv, it);
    01770                 }
    01771                 if (it != it_end && (*it).first == layout_type::index_m (i, j))
    01772                     return const_iterator2 (*this, rank, i, j, itv, it);
    01773                 if (direction > 0) {
    01774                     if (layout_type::fast_j ()) {
    01775                         if (it == it_end)
    01776                             return const_iterator2 (*this, rank, i, j, itv, it);
    01777                         j = (*it).first;
    01778                     } else {
    01779                         if (j >= size2_)
    01780                             return const_iterator2 (*this, rank, i, j, itv, it);
    01781                         ++ j;
    01782                     }
    01783                 } else /* if (direction < 0)  */ {
    01784                     if (layout_type::fast_j ()) {
    01785                         if (it == (*itv).second.begin ())
    01786                             return const_iterator2 (*this, rank, i, j, itv, it);
    01787                         -- it;
    01788                         j = (*it).first;
    01789                     } else {
    01790                         if (j == 0)
    01791                             return const_iterator2 (*this, rank, i, j, itv, it);
    01792                         -- j;
    01793                     }
    01794                 }
    01795             }
    01796         }
    01797         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    01798         iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) {
    01799             BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ());
    01800             for (;;) {
    01801                 vector_subiterator_type itv (data ().lower_bound (layout_type::index_M (i, j)));
    01802                 vector_subiterator_type itv_end (data ().end ());
    01803                 if (itv == itv_end)
    01804                     return iterator2 (*this, rank, i, j, itv_end, (*(-- itv)).second.end ());
    01805
    01806                 subiterator_type it ((*itv).second.lower_bound (layout_type::index_m (i, j)));
    01807                 subiterator_type it_end ((*itv).second.end ());
    01808                 if (rank == 0) {
    01809                     // advance to the first available major index
    01810                     size_type M = itv->first;
    01811                     size_type m;
    01812                     if (it != it_end) {
    01813                         m = it->first;
    01814                     } else {
    01815                         m = layout_type::size_m(size1_, size2_);
    01816                     }
    01817                     size_type first_j = layout_type::index_m(M,m);
    01818                     return iterator2 (*this, rank, i, first_j, itv, it);
    01819                 }
    01820                 if (it != it_end && (*it).first == layout_type::index_m (i, j))
    01821                     return iterator2 (*this, rank, i, j, itv, it);
    01822                 if (direction > 0) {
    01823                     if (layout_type::fast_j ()) {
    01824                         if (it == it_end)
    01825                             return iterator2 (*this, rank, i, j, itv, it);
    01826                         j = (*it).first;
    01827                     } else {
    01828                         if (j >= size2_)
    01829                             return iterator2 (*this, rank, i, j, itv, it);
    01830                         ++ j;
    01831                     }
    01832                 } else /* if (direction < 0)  */ {
    01833                     if (layout_type::fast_j ()) {
    01834                         if (it == (*itv).second.begin ())
    01835                             return iterator2 (*this, rank, i, j, itv, it);
    01836                         -- it;
    01837                         j = (*it).first;
    01838                     } else {
    01839                         if (j == 0)
    01840                             return iterator2 (*this, rank, i, j, itv, it);
    01841                         -- j;
    01842                     }
    01843                 }
    01844             }
    01845         }
    01846
    01847         class const_iterator1:
    01848             public container_const_reference<mapped_vector_of_mapped_vector>,
    01849             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
    01850                                                const_iterator1, value_type> {
    01851         public:
    01852             typedef typename mapped_vector_of_mapped_vector::value_type value_type;
    01853             typedef typename mapped_vector_of_mapped_vector::difference_type difference_type;
    01854             typedef typename mapped_vector_of_mapped_vector::const_reference reference;
    01855             typedef const typename mapped_vector_of_mapped_vector::pointer pointer;
    01856
    01857             typedef const_iterator2 dual_iterator_type;
    01858             typedef const_reverse_iterator2 dual_reverse_iterator_type;
    01859
    01860             // Construction and destruction
    01861             BOOST_UBLAS_INLINE
    01862             const_iterator1 ():
    01863                 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
    01864             BOOST_UBLAS_INLINE
    01865             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):
    01866                 container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
    01867             BOOST_UBLAS_INLINE
    01868             const_iterator1 (const iterator1 &it):
    01869                 container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {}
    01870
    01871             // Arithmetic
    01872             BOOST_UBLAS_INLINE
    01873             const_iterator1 &operator ++ () {
    01874                 if (rank_ == 1 && layout_type::fast_i ())
    01875                     ++ it_;
    01876                 else {
    01877                     const self_type &m = (*this) ();
    01878                     if (rank_ == 0) {
    01879                         ++ itv_;
    01880                         i_ = itv_->first;
    01881                     } else {
    01882                         i_ = index1 () + 1;
    01883                     }
    01884                     if (rank_ == 1 && ++ itv_ == m.end1 ().itv_)
    01885                         *this = m.find1 (rank_, i_, j_, 1);
    01886                     else if (rank_ == 1) {
    01887                         it_ = (*itv_).second.begin ();
    01888                         if (it_ == (*itv_).second.end () || index2 () != j_)
    01889                             *this = m.find1 (rank_, i_, j_, 1);
    01890                     }
    01891                 }
    01892                 return *this;
    01893             }
    01894             BOOST_UBLAS_INLINE
    01895             const_iterator1 &operator -- () {
    01896                 if (rank_ == 1 && layout_type::fast_i ())
    01897                     -- it_;
    01898                 else {
    01899                     const self_type &m = (*this) ();
    01900                     if (rank_ == 0) {
    01901                         -- itv_;
    01902                         i_ = itv_->first;
    01903                     } else {
    01904                         i_ = index1 () - 1;
    01905                     }
    01906                     // FIXME: this expression should never become true!
    01907                     if (rank_ == 1 && -- itv_ == m.end1 ().itv_)
    01908                         *this = m.find1 (rank_, i_, j_, -1);
    01909                     else if (rank_ == 1) {
    01910                         it_ = (*itv_).second.begin ();
    01911                         if (it_ == (*itv_).second.end () || index2 () != j_)
    01912                             *this = m.find1 (rank_, i_, j_, -1);
    01913                     }
    01914                 }
    01915                 return *this;
    01916             }
    01917
    01918             // Dereference
    01919             BOOST_UBLAS_INLINE
    01920             const_reference operator * () const {
    01921                 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
    01922                 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
    01923                 if (rank_ == 1) {
    01924                     return (*it_).second;
    01925                 } else {
    01926                     return (*this) () (i_, j_);
    01927                 }
    01928             }
    01929
    01930 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    01931             BOOST_UBLAS_INLINE
    01932 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    01933             typename self_type::
    01934 #endif
    01935             const_iterator2 begin () const {
    01936                 const self_type &m = (*this) ();
    01937                 return m.find2 (1, index1 (), 0);
    01938             }
    01939             BOOST_UBLAS_INLINE
    01940 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    01941             typename self_type::
    01942 #endif
    01943             const_iterator2 end () const {
    01944                 const self_type &m = (*this) ();
    01945                 return m.find2 (1, index1 (), m.size2 ());
    01946             }
    01947             BOOST_UBLAS_INLINE
    01948 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    01949             typename self_type::
    01950 #endif
    01951             const_reverse_iterator2 rbegin () const {
    01952                 return const_reverse_iterator2 (end ());
    01953             }
    01954             BOOST_UBLAS_INLINE
    01955 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    01956             typename self_type::
    01957 #endif
    01958             const_reverse_iterator2 rend () const {
    01959                 return const_reverse_iterator2 (begin ());
    01960             }
    01961 #endif
    01962
    01963             // Indices
    01964             BOOST_UBLAS_INLINE
    01965             size_type index1 () const {
    01966                 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
    01967                 if (rank_ == 1) {
    01968                     BOOST_UBLAS_CHECK (layout_type::index_M ((*itv_).first, (*it_).first) < (*this) ().size1 (), bad_index ());
    01969                     return layout_type::index_M ((*itv_).first, (*it_).first);
    01970                 } else {
    01971                     return i_;
    01972                 }
    01973             }
    01974             BOOST_UBLAS_INLINE
    01975             size_type index2 () const {
    01976                 if (rank_ == 1) {
    01977                     BOOST_UBLAS_CHECK (layout_type::index_m ((*itv_).first, (*it_).first) < (*this) ().size2 (), bad_index ());
    01978                     return layout_type::index_m ((*itv_).first, (*it_).first);
    01979                 } else {
    01980                     return j_;
    01981                 }
    01982             }
    01983
    01984             // Assignment
    01985             BOOST_UBLAS_INLINE
    01986             const_iterator1 &operator = (const const_iterator1 &it) {
    01987                 container_const_reference<self_type>::assign (&it ());
    01988                 rank_ = it.rank_;
    01989                 i_ = it.i_;
    01990                 j_ = it.j_;
    01991                 itv_ = it.itv_;
    01992                 it_ = it.it_;
    01993                 return *this;
    01994             }
    01995
    01996             // Comparison
    01997             BOOST_UBLAS_INLINE
    01998             bool operator == (const const_iterator1 &it) const {
    01999                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
    02000                 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
    02001                 if (rank_ == 1 || it.rank_ == 1) {
    02002                     return it_ == it.it_;
    02003                 } else {
    02004                     return i_ == it.i_ && j_ == it.j_;
    02005                 }
    02006             }
    02007
    02008         private:
    02009             int rank_;
    02010             size_type i_;
    02011             size_type j_;
    02012             vector_const_subiterator_type itv_;
    02013             const_subiterator_type it_;
    02014         };
    02015
    02016         BOOST_UBLAS_INLINE
    02017         const_iterator1 begin1 () const {
    02018             return find1 (0, 0, 0);
    02019         }
    02020         BOOST_UBLAS_INLINE
    02021         const_iterator1 end1 () const {
    02022             return find1 (0, size1_, 0);
    02023         }
    02024
    02025         class iterator1:
    02026             public container_reference<mapped_vector_of_mapped_vector>,
    02027             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
    02028                                                iterator1, value_type> {
    02029         public:
    02030             typedef typename mapped_vector_of_mapped_vector::value_type value_type;
    02031             typedef typename mapped_vector_of_mapped_vector::difference_type difference_type;
    02032             typedef typename mapped_vector_of_mapped_vector::true_reference reference;
    02033             typedef typename mapped_vector_of_mapped_vector::pointer pointer;
    02034
    02035             typedef iterator2 dual_iterator_type;
    02036             typedef reverse_iterator2 dual_reverse_iterator_type;
    02037
    02038             // Construction and destruction
    02039             BOOST_UBLAS_INLINE
    02040             iterator1 ():
    02041                 container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
    02042             BOOST_UBLAS_INLINE
    02043             iterator1 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it):
    02044                 container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
    02045
    02046             // Arithmetic
    02047             BOOST_UBLAS_INLINE
    02048             iterator1 &operator ++ () {
    02049                 if (rank_ == 1 && layout_type::fast_i ())
    02050                     ++ it_;
    02051                 else {
    02052                     self_type &m = (*this) ();
    02053                     if (rank_ == 0) {
    02054                         ++ itv_;
    02055                         i_ = itv_->first;
    02056                     } else {
    02057                         i_ = index1 () + 1;
    02058                     }
    02059                     if (rank_ == 1 && ++ itv_ == m.end1 ().itv_)
    02060                         *this = m.find1 (rank_, i_, j_, 1);
    02061                     else if (rank_ == 1) {
    02062                         it_ = (*itv_).second.begin ();
    02063                         if (it_ == (*itv_).second.end () || index2 () != j_)
    02064                             *this = m.find1 (rank_, i_, j_, 1);
    02065                     }
    02066                 }
    02067                 return *this;
    02068             }
    02069             BOOST_UBLAS_INLINE
    02070             iterator1 &operator -- () {
    02071                 if (rank_ == 1 && layout_type::fast_i ())
    02072                     -- it_;
    02073                 else {
    02074                     self_type &m = (*this) ();
    02075                     if (rank_ == 0) {
    02076                         -- itv_;
    02077                         i_ = itv_->first;
    02078                     } else {
    02079                         i_ = index1 () - 1;
    02080                     }
    02081                     // FIXME: this expression should never become true!
    02082                     if (rank_ == 1 && -- itv_ == m.end1 ().itv_)
    02083                         *this = m.find1 (rank_, i_, j_, -1);
    02084                     else if (rank_ == 1) {
    02085                         it_ = (*itv_).second.begin ();
    02086                         if (it_ == (*itv_).second.end () || index2 () != j_)
    02087                             *this = m.find1 (rank_, i_, j_, -1);
    02088                     }
    02089                 }
    02090                 return *this;
    02091             }
    02092
    02093             // Dereference
    02094             BOOST_UBLAS_INLINE
    02095             reference operator * () const {
    02096                 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
    02097                 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
    02098                 if (rank_ == 1) {
    02099                     return (*it_).second;
    02100                 } else {
    02101                     return (*this) ().at_element (i_, j_);
    02102                 }
    02103             }
    02104
    02105 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    02106             BOOST_UBLAS_INLINE
    02107 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    02108             typename self_type::
    02109 #endif
    02110             iterator2 begin () const {
    02111                 self_type &m = (*this) ();
    02112                 return m.find2 (1, index1 (), 0);
    02113             }
    02114             BOOST_UBLAS_INLINE
    02115 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    02116             typename self_type::
    02117 #endif
    02118             iterator2 end () const {
    02119                 self_type &m = (*this) ();
    02120                 return m.find2 (1, index1 (), m.size2 ());
    02121             }
    02122             BOOST_UBLAS_INLINE
    02123 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    02124             typename self_type::
    02125 #endif
    02126             reverse_iterator2 rbegin () const {
    02127                 return reverse_iterator2 (end ());
    02128             }
    02129             BOOST_UBLAS_INLINE
    02130 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    02131             typename self_type::
    02132 #endif
    02133             reverse_iterator2 rend () const {
    02134                 return reverse_iterator2 (begin ());
    02135             }
    02136 #endif
    02137
    02138             // Indices
    02139             BOOST_UBLAS_INLINE
    02140             size_type index1 () const {
    02141                 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
    02142                 if (rank_ == 1) {
    02143                     BOOST_UBLAS_CHECK (layout_type::index_M ((*itv_).first, (*it_).first) < (*this) ().size1 (), bad_index ());
    02144                     return layout_type::index_M ((*itv_).first, (*it_).first);
    02145                 } else {
    02146                     return i_;
    02147                 }
    02148             }
    02149             BOOST_UBLAS_INLINE
    02150             size_type index2 () const {
    02151                 if (rank_ == 1) {
    02152                     BOOST_UBLAS_CHECK (layout_type::index_m ((*itv_).first, (*it_).first) < (*this) ().size2 (), bad_index ());
    02153                     return layout_type::index_m ((*itv_).first, (*it_).first);
    02154                 } else {
    02155                     return j_;
    02156                 }
    02157             }
    02158
    02159             // Assignment
    02160             BOOST_UBLAS_INLINE
    02161             iterator1 &operator = (const iterator1 &it) {
    02162                 container_reference<self_type>::assign (&it ());
    02163                 rank_ = it.rank_;
    02164                 i_ = it.i_;
    02165                 j_ = it.j_;
    02166                 itv_ = it.itv_;
    02167                 it_ = it.it_;
    02168                 return *this;
    02169             }
    02170
    02171             // Comparison
    02172             BOOST_UBLAS_INLINE
    02173             bool operator == (const iterator1 &it) const {
    02174                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
    02175                 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
    02176                 if (rank_ == 1 || it.rank_ == 1) {
    02177                     return it_ == it.it_;
    02178                 } else {
    02179                     return i_ == it.i_ && j_ == it.j_;
    02180                 }
    02181             }
    02182
    02183         private:
    02184             int rank_;
    02185             size_type i_;
    02186             size_type j_;
    02187             vector_subiterator_type itv_;
    02188             subiterator_type it_;
    02189
    02190             friend class const_iterator1;
    02191         };
    02192
    02193         BOOST_UBLAS_INLINE
    02194         iterator1 begin1 () {
    02195             return find1 (0, 0, 0);
    02196         }
    02197         BOOST_UBLAS_INLINE
    02198         iterator1 end1 () {
    02199             return find1 (0, size1_, 0);
    02200         }
    02201
    02202         class const_iterator2:
    02203             public container_const_reference<mapped_vector_of_mapped_vector>,
    02204             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
    02205                                                const_iterator2, value_type> {
    02206         public:
    02207             typedef typename mapped_vector_of_mapped_vector::value_type value_type;
    02208             typedef typename mapped_vector_of_mapped_vector::difference_type difference_type;
    02209             typedef typename mapped_vector_of_mapped_vector::const_reference reference;
    02210             typedef const typename mapped_vector_of_mapped_vector::pointer pointer;
    02211
    02212             typedef const_iterator1 dual_iterator_type;
    02213             typedef const_reverse_iterator1 dual_reverse_iterator_type;
    02214
    02215             // Construction and destruction
    02216             BOOST_UBLAS_INLINE
    02217             const_iterator2 ():
    02218                 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
    02219             BOOST_UBLAS_INLINE
    02220             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):
    02221                 container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
    02222             BOOST_UBLAS_INLINE
    02223             const_iterator2 (const iterator2 &it):
    02224                 container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {}
    02225
    02226             // Arithmetic
    02227             BOOST_UBLAS_INLINE
    02228             const_iterator2 &operator ++ () {
    02229                 if (rank_ == 1 && layout_type::fast_j ())
    02230                     ++ it_;
    02231                 else {
    02232                     const self_type &m = (*this) ();
    02233                     if (rank_ == 0) {
    02234                         ++ itv_;
    02235                         j_ = itv_->first;
    02236                     } else {
    02237                         j_ = index2 () + 1;
    02238                     }
    02239                     if (rank_ == 1 && ++ itv_ == m.end2 ().itv_)
    02240                         *this = m.find2 (rank_, i_, j_, 1);
    02241                     else if (rank_ == 1) {
    02242                         it_ = (*itv_).second.begin ();
    02243                         if (it_ == (*itv_).second.end () || index1 () != i_)
    02244                             *this = m.find2 (rank_, i_, j_, 1);
    02245                     }
    02246                 }
    02247                 return *this;
    02248             }
    02249             BOOST_UBLAS_INLINE
    02250             const_iterator2 &operator -- () {
    02251                 if (rank_ == 1 && layout_type::fast_j ())
    02252                     -- it_;
    02253                 else {
    02254                     const self_type &m = (*this) ();
    02255                     if (rank_ == 0) {
    02256                         -- itv_;
    02257                         j_ = itv_->first;
    02258                     } else {
    02259                         j_ = index2 () - 1;
    02260                     }
    02261                     // FIXME: this expression should never become true!
    02262                     if (rank_ == 1 && -- itv_ == m.end2 ().itv_)
    02263                         *this = m.find2 (rank_, i_, j_, -1);
    02264                     else if (rank_ == 1) {
    02265                         it_ = (*itv_).second.begin ();
    02266                         if (it_ == (*itv_).second.end () || index1 () != i_)
    02267                             *this = m.find2 (rank_, i_, j_, -1);
    02268                     }
    02269                 }
    02270                 return *this;
    02271             }
    02272
    02273             // Dereference
    02274             BOOST_UBLAS_INLINE
    02275             const_reference operator * () const {
    02276                 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
    02277                 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
    02278                 if (rank_ == 1) {
    02279                     return (*it_).second;
    02280                 } else {
    02281                     return (*this) () (i_, j_);
    02282                 }
    02283             }
    02284
    02285 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    02286             BOOST_UBLAS_INLINE
    02287 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    02288             typename self_type::
    02289 #endif
    02290             const_iterator1 begin () const {
    02291                 const self_type &m = (*this) ();
    02292                 return m.find1 (1, 0, index2 ());
    02293             }
    02294             BOOST_UBLAS_INLINE
    02295 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    02296             typename self_type::
    02297 #endif
    02298             const_iterator1 end () const {
    02299                 const self_type &m = (*this) ();
    02300                 return m.find1 (1, m.size1 (), index2 ());
    02301             }
    02302             BOOST_UBLAS_INLINE
    02303 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    02304             typename self_type::
    02305 #endif
    02306             const_reverse_iterator1 rbegin () const {
    02307                 return const_reverse_iterator1 (end ());
    02308             }
    02309             BOOST_UBLAS_INLINE
    02310 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    02311             typename self_type::
    02312 #endif
    02313             const_reverse_iterator1 rend () const {
    02314                 return const_reverse_iterator1 (begin ());
    02315             }
    02316 #endif
    02317
    02318             // Indices
    02319             BOOST_UBLAS_INLINE
    02320             size_type index1 () const {
    02321                 if (rank_ == 1) {
    02322                     BOOST_UBLAS_CHECK (layout_type::index_M ((*itv_).first, (*it_).first) < (*this) ().size1 (), bad_index ());
    02323                     return layout_type::index_M ((*itv_).first, (*it_).first);
    02324                 } else {
    02325                     return i_;
    02326                 }
    02327             }
    02328             BOOST_UBLAS_INLINE
    02329             size_type index2 () const {
    02330                 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
    02331                 if (rank_ == 1) {
    02332                     BOOST_UBLAS_CHECK (layout_type::index_m ((*itv_).first, (*it_).first) < (*this) ().size2 (), bad_index ());
    02333                     return layout_type::index_m ((*itv_).first, (*it_).first);
    02334                 } else {
    02335                     return j_;
    02336                 }
    02337             }
    02338
    02339             // Assignment
    02340             BOOST_UBLAS_INLINE
    02341             const_iterator2 &operator = (const const_iterator2 &it) {
    02342                 container_const_reference<self_type>::assign (&it ());
    02343                 rank_ = it.rank_;
    02344                 i_ = it.i_;
    02345                 j_ = it.j_;
    02346                 itv_ = it.itv_;
    02347                 it_ = it.it_;
    02348                 return *this;
    02349             }
    02350
    02351             // Comparison
    02352             BOOST_UBLAS_INLINE
    02353             bool operator == (const const_iterator2 &it) const {
    02354                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
    02355                 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
    02356                 if (rank_ == 1 || it.rank_ == 1) {
    02357                     return it_ == it.it_;
    02358                 } else {
    02359                     return i_ == it.i_ && j_ == it.j_;
    02360                 }
    02361             }
    02362
    02363         private:
    02364             int rank_;
    02365             size_type i_;
    02366             size_type j_;
    02367             vector_const_subiterator_type itv_;
    02368             const_subiterator_type it_;
    02369         };
    02370
    02371         BOOST_UBLAS_INLINE
    02372         const_iterator2 begin2 () const {
    02373             return find2 (0, 0, 0);
    02374         }
    02375         BOOST_UBLAS_INLINE
    02376         const_iterator2 end2 () const {
    02377             return find2 (0, 0, size2_);
    02378         }
    02379
    02380         class iterator2:
    02381             public container_reference<mapped_vector_of_mapped_vector>,
    02382             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
    02383                                                iterator2, value_type> {
    02384         public:
    02385             typedef typename mapped_vector_of_mapped_vector::value_type value_type;
    02386             typedef typename mapped_vector_of_mapped_vector::difference_type difference_type;
    02387             typedef typename mapped_vector_of_mapped_vector::true_reference reference;
    02388             typedef typename mapped_vector_of_mapped_vector::pointer pointer;
    02389
    02390             typedef iterator1 dual_iterator_type;
    02391             typedef reverse_iterator1 dual_reverse_iterator_type;
    02392
    02393             // Construction and destruction
    02394             BOOST_UBLAS_INLINE
    02395             iterator2 ():
    02396                 container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
    02397             BOOST_UBLAS_INLINE
    02398             iterator2 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it):
    02399                 container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
    02400
    02401             // Arithmetic
    02402             BOOST_UBLAS_INLINE
    02403             iterator2 &operator ++ () {
    02404                 if (rank_ == 1 && layout_type::fast_j ())
    02405                     ++ it_;
    02406                 else {
    02407                     self_type &m = (*this) ();
    02408                     if (rank_ == 0) {
    02409                         ++ itv_;
    02410                         j_ = itv_->first;
    02411                     } else {
    02412                         j_ = index2 () + 1;
    02413                     }
    02414                     if (rank_ == 1 && ++ itv_ == m.end2 ().itv_)
    02415                         *this = m.find2 (rank_, i_, j_, 1);
    02416                     else if (rank_ == 1) {
    02417                         it_ = (*itv_).second.begin ();
    02418                         if (it_ == (*itv_).second.end () || index1 () != i_)
    02419                             *this = m.find2 (rank_, i_, j_, 1);
    02420                     }
    02421                 }
    02422                 return *this;
    02423             }
    02424             BOOST_UBLAS_INLINE
    02425             iterator2 &operator -- () {
    02426                 if (rank_ == 1 && layout_type::fast_j ())
    02427                     -- it_;
    02428                 else {
    02429                     self_type &m = (*this) ();
    02430                     if (rank_ == 0) {
    02431                         -- itv_;
    02432                         j_ = itv_->first;
    02433                     } else {
    02434                         j_ = index2 () - 1;
    02435                     }
    02436                     // FIXME: this expression should never become true!
    02437                     if (rank_ == 1 && -- itv_ == m.end2 ().itv_)
    02438                         *this = m.find2 (rank_, i_, j_, -1);
    02439                     else if (rank_ == 1) {
    02440                         it_ = (*itv_).second.begin ();
    02441                         if (it_ == (*itv_).second.end () || index1 () != i_)
    02442                             *this = m.find2 (rank_, i_, j_, -1);
    02443                     }
    02444                 }
    02445                 return *this;
    02446             }
    02447
    02448             // Dereference
    02449             BOOST_UBLAS_INLINE
    02450             reference operator * () const {
    02451                 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
    02452                 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
    02453                 if (rank_ == 1) {
    02454                     return (*it_).second;
    02455                 } else {
    02456                     return (*this) ().at_element (i_, j_);
    02457                 }
    02458             }
    02459
    02460 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    02461             BOOST_UBLAS_INLINE
    02462 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    02463             typename self_type::
    02464 #endif
    02465             iterator1 begin () const {
    02466                 self_type &m = (*this) ();
    02467                 return m.find1 (1, 0, index2 ());
    02468             }
    02469             BOOST_UBLAS_INLINE
    02470 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    02471             typename self_type::
    02472 #endif
    02473             iterator1 end () const {
    02474                 self_type &m = (*this) ();
    02475                 return m.find1 (1, m.size1 (), index2 ());
    02476             }
    02477             BOOST_UBLAS_INLINE
    02478 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    02479             typename self_type::
    02480 #endif
    02481             reverse_iterator1 rbegin () const {
    02482                 return reverse_iterator1 (end ());
    02483             }
    02484             BOOST_UBLAS_INLINE
    02485 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    02486             typename self_type::
    02487 #endif
    02488             reverse_iterator1 rend () const {
    02489                 return reverse_iterator1 (begin ());
    02490             }
    02491 #endif
    02492
    02493             // Indices
    02494             BOOST_UBLAS_INLINE
    02495             size_type index1 () const {
    02496                 if (rank_ == 1) {
    02497                     BOOST_UBLAS_CHECK (layout_type::index_M ((*itv_).first, (*it_).first) < (*this) ().size1 (), bad_index ());
    02498                     return layout_type::index_M ((*itv_).first, (*it_).first);
    02499                 } else {
    02500                     return i_;
    02501                 }
    02502             }
    02503             BOOST_UBLAS_INLINE
    02504             size_type index2 () const {
    02505                 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
    02506                 if (rank_ == 1) {
    02507                     BOOST_UBLAS_CHECK (layout_type::index_m ((*itv_).first, (*it_).first) < (*this) ().size2 (), bad_index ());
    02508                     return layout_type::index_m ((*itv_).first, (*it_).first);
    02509                 } else {
    02510                     return j_;
    02511                 }
    02512             }
    02513
    02514             // Assignment
    02515             BOOST_UBLAS_INLINE
    02516             iterator2 &operator = (const iterator2 &it) {
    02517                 container_reference<self_type>::assign (&it ());
    02518                 rank_ = it.rank_;
    02519                 i_ = it.i_;
    02520                 j_ = it.j_;
    02521                 itv_ = it.itv_;
    02522                 it_ = it.it_;
    02523                 return *this;
    02524             }
    02525
    02526             // Comparison
    02527             BOOST_UBLAS_INLINE
    02528             bool operator == (const iterator2 &it) const {
    02529                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
    02530                 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
    02531                 if (rank_ == 1 || it.rank_ == 1) {
    02532                     return it_ == it.it_;
    02533                 } else {
    02534                     return i_ == it.i_ && j_ == it.j_;
    02535                 }
    02536             }
    02537
    02538         private:
    02539             int rank_;
    02540             size_type i_;
    02541             size_type j_;
    02542             vector_subiterator_type itv_;
    02543             subiterator_type it_;
    02544
    02545             friend class const_iterator2;
    02546         };
    02547
    02548         BOOST_UBLAS_INLINE
    02549         iterator2 begin2 () {
    02550             return find2 (0, 0, 0);
    02551         }
    02552         BOOST_UBLAS_INLINE
    02553         iterator2 end2 () {
    02554             return find2 (0, 0, size2_);
    02555         }
    02556
    02557         // Reverse iterators
    02558
    02559         BOOST_UBLAS_INLINE
    02560         const_reverse_iterator1 rbegin1 () const {
    02561             return const_reverse_iterator1 (end1 ());
    02562         }
    02563         BOOST_UBLAS_INLINE
    02564         const_reverse_iterator1 rend1 () const {
    02565             return const_reverse_iterator1 (begin1 ());
    02566         }
    02567
    02568         BOOST_UBLAS_INLINE
    02569         reverse_iterator1 rbegin1 () {
    02570             return reverse_iterator1 (end1 ());
    02571         }
    02572         BOOST_UBLAS_INLINE
    02573         reverse_iterator1 rend1 () {
    02574             return reverse_iterator1 (begin1 ());
    02575         }
    02576
    02577         BOOST_UBLAS_INLINE
    02578         const_reverse_iterator2 rbegin2 () const {
    02579             return const_reverse_iterator2 (end2 ());
    02580         }
    02581         BOOST_UBLAS_INLINE
    02582         const_reverse_iterator2 rend2 () const {
    02583             return const_reverse_iterator2 (begin2 ());
    02584         }
    02585
    02586         BOOST_UBLAS_INLINE
    02587         reverse_iterator2 rbegin2 () {
    02588             return reverse_iterator2 (end2 ());
    02589         }
    02590         BOOST_UBLAS_INLINE
    02591         reverse_iterator2 rend2 () {
    02592             return reverse_iterator2 (begin2 ());
    02593         }
    02594
    02595          // Serialization
    02596         template<class Archive>
    02597         void serialize(Archive & ar, const unsigned int /* file_version */){
    02598             serialization::collection_size_type s1 (size1_);
    02599             serialization::collection_size_type s2 (size2_);
    02600             ar & serialization::make_nvp("size1",s1);
    02601             ar & serialization::make_nvp("size2",s2);
    02602             if (Archive::is_loading::value) {
    02603                 size1_ = s1;
    02604                 size2_ = s2;
    02605             }
    02606             ar & serialization::make_nvp("data", data_);
    02607         }
    02608
    02609     private:
    02610         size_type size1_;
    02611         size_type size2_;
    02612         array_type data_;
    02613         static const value_type zero_;
    02614     };
    02615
    02616     template<class T, class L, class A>
    02617     const typename mapped_vector_of_mapped_vector<T, L, A>::value_type mapped_vector_of_mapped_vector<T, L, A>::zero_ = value_type/*zero*/();
    02618
    02619
    02620     // Comperssed array based sparse matrix class
    02621     // Thanks to Kresimir Fresl for extending this to cover different index bases.
    02622     template<class T, class L, std::size_t IB, class IA, class TA>
    02623     class compressed_matrix:
    02624         public matrix_container<compressed_matrix<T, L, IB, IA, TA> > {
    02625
    02626         typedef T &true_reference;
    02627         typedef T *pointer;
    02628         typedef const T *const_pointer;
    02629         typedef L layout_type;
    02630         typedef compressed_matrix<T, L, IB, IA, TA> self_type;
    02631     public:
    02632 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
    02633         using matrix_container<self_type>::operator ();
    02634 #endif
    02635         // ISSUE require type consistency check
    02636         // is_convertable (IA::size_type, TA::size_type)
    02637         typedef typename IA::value_type size_type;
    02638         // size_type for the data arrays.
    02639         typedef typename IA::size_type array_size_type;
    02640         // FIXME difference type for sparse storage iterators should it be in the container?
    02641         typedef typename IA::difference_type difference_type;
    02642         typedef T value_type;
    02643         typedef const T &const_reference;
    02644 #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
    02645         typedef T &reference;
    02646 #else
    02647         typedef sparse_matrix_element<self_type> reference;
    02648 #endif
    02649         typedef IA index_array_type;
    02650         typedef TA value_array_type;
    02651         typedef const matrix_reference<const self_type> const_closure_type;
    02652         typedef matrix_reference<self_type> closure_type;
    02653         typedef compressed_vector<T, IB, IA, TA> vector_temporary_type;
    02654         typedef self_type matrix_temporary_type;
    02655         typedef sparse_tag storage_category;
    02656         typedef typename L::orientation_category orientation_category;
    02657
    02658         // Construction and destruction
    02659         BOOST_UBLAS_INLINE
    02660         compressed_matrix ():
    02661             matrix_container<self_type> (),
    02662             size1_ (0), size2_ (0), capacity_ (restrict_capacity (0)),
    02663             filled1_ (1), filled2_ (0),
    02664             index1_data_ (layout_type::size_M (size1_, size2_) + 1), index2_data_ (capacity_), value_data_ (capacity_) {
    02665             index1_data_ [filled1_ - 1] = k_based (filled2_);
    02666             storage_invariants ();
    02667         }
    02668         BOOST_UBLAS_INLINE
    02669         compressed_matrix (size_type size1, size_type size2, size_type non_zeros = 0):
    02670             matrix_container<self_type> (),
    02671             size1_ (size1), size2_ (size2), capacity_ (restrict_capacity (non_zeros)),
    02672             filled1_ (1), filled2_ (0),
    02673             index1_data_ (layout_type::size_M (size1_, size2_) + 1), index2_data_ (capacity_), value_data_ (capacity_) {
    02674             index1_data_ [filled1_ - 1] = k_based (filled2_);
    02675             storage_invariants ();
    02676         }
    02677         BOOST_UBLAS_INLINE
    02678         compressed_matrix (const compressed_matrix &m):
    02679             matrix_container<self_type> (),
    02680             size1_ (m.size1_), size2_ (m.size2_), capacity_ (m.capacity_),
    02681             filled1_ (m.filled1_), filled2_ (m.filled2_),
    02682             index1_data_ (m.index1_data_), index2_data_ (m.index2_data_), value_data_ (m.value_data_) {
    02683             storage_invariants ();
    02684         }
    02685
    02686         BOOST_UBLAS_INLINE
    02687         compressed_matrix (const coordinate_matrix<T, L, IB, IA, TA> &m):
    02688             matrix_container<self_type> (),
    02689             size1_ (m.size1()), size2_ (m.size2()),
    02690             index1_data_ (layout_type::size_M (size1_, size2_) + 1)
    02691         {
    02692             m.sort();
    02693             reserve(m.nnz(), false);
    02694             filled2_ = m.nnz();
    02695             const_subiterator_type  i_start = m.index1_data().begin();
    02696             const_subiterator_type  i_end   = (i_start + filled2_);
    02697             const_subiterator_type  i = i_start;
    02698             size_type r = 1;
    02699             for (; (r < layout_type::size_M (size1_, size2_)) && (i != i_end); ++r) {
    02700                 i = std::lower_bound(i, i_end, r);
    02701                 index1_data_[r] = k_based( i - i_start );
    02702             }
    02703             filled1_ = r + 1;
    02704             std::copy( m.index2_data().begin(), m.index2_data().begin() + filled2_, index2_data_.begin());
    02705             std::copy( m.value_data().begin(), m.value_data().begin() + filled2_, value_data_.begin());
    02706             index1_data_ [filled1_ - 1] = k_based(filled2_);
    02707             storage_invariants ();
    02708         }
    02709
    02710        template<class AE>
    02711        BOOST_UBLAS_INLINE
    02712        compressed_matrix (const matrix_expression<AE> &ae, size_type non_zeros = 0):
    02713             matrix_container<self_type> (),
    02714             size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), capacity_ (restrict_capacity (non_zeros)),
    02715             filled1_ (1), filled2_ (0),
    02716             index1_data_ (layout_type::size_M (ae ().size1 (), ae ().size2 ()) + 1),
    02717             index2_data_ (capacity_), value_data_ (capacity_) {
    02718             index1_data_ [filled1_ - 1] = k_based (filled2_);
    02719             storage_invariants ();
    02720             matrix_assign<scalar_assign> (*this, ae);
    02721         }
    02722
    02723         // Accessors
    02724         BOOST_UBLAS_INLINE
    02725         size_type size1 () const {
    02726             return size1_;
    02727         }
    02728         BOOST_UBLAS_INLINE
    02729         size_type size2 () const {
    02730             return size2_;
    02731         }
    02732         BOOST_UBLAS_INLINE
    02733         size_type nnz_capacity () const {
    02734             return capacity_;
    02735         }
    02736         BOOST_UBLAS_INLINE
    02737         size_type nnz () const {
    02738             return filled2_;
    02739         }
    02740
    02741         // Storage accessors
    02742         BOOST_UBLAS_INLINE
    02743         static size_type index_base () {
    02744             return IB;
    02745         }
    02746         BOOST_UBLAS_INLINE
    02747         array_size_type filled1 () const {
    02748             return filled1_;
    02749         }
    02750         BOOST_UBLAS_INLINE
    02751         array_size_type filled2 () const {
    02752             return filled2_;
    02753         }
    02754         BOOST_UBLAS_INLINE
    02755         const index_array_type &index1_data () const {
    02756             return index1_data_;
    02757         }
    02758         BOOST_UBLAS_INLINE
    02759         const index_array_type &index2_data () const {
    02760             return index2_data_;
    02761         }
    02762         BOOST_UBLAS_INLINE
    02763         const value_array_type &value_data () const {
    02764             return value_data_;
    02765         }
    02766         BOOST_UBLAS_INLINE
    02767         void set_filled (const array_size_type& filled1, const array_size_type& filled2) {
    02768             filled1_ = filled1;
    02769             filled2_ = filled2;
    02770             storage_invariants ();
    02771         }
    02772         BOOST_UBLAS_INLINE
    02773         index_array_type &index1_data () {
    02774             return index1_data_;
    02775         }
    02776         BOOST_UBLAS_INLINE
    02777         index_array_type &index2_data () {
    02778             return index2_data_;
    02779         }
    02780         BOOST_UBLAS_INLINE
    02781         value_array_type &value_data () {
    02782             return value_data_;
    02783         }
    02784         BOOST_UBLAS_INLINE
    02785         void complete_index1_data () {
    02786             while (filled1_ <= layout_type::size_M (size1_, size2_)) {
    02787                 this->index1_data_ [filled1_] = k_based (filled2_);
    02788                 ++ this->filled1_;
    02789             }
    02790         }
    02791
    02792         // Resizing
    02793     private:
    02794         BOOST_UBLAS_INLINE
    02795         size_type restrict_capacity (size_type non_zeros) const {
    02796             non_zeros = (std::max) (non_zeros, (std::min) (size1_, size2_));
    02797             // Guarding against overflow - Thanks to Alexei Novakov for the hint.
    02798             // non_zeros = (std::min) (non_zeros, size1_ * size2_);
    02799             if (size1_ > 0 && non_zeros / size1_ >= size2_)
    02800                 non_zeros = size1_ * size2_;
    02801             return non_zeros;
    02802         }
    02803     public:
    02804         BOOST_UBLAS_INLINE
    02805         void resize (size_type size1, size_type size2, bool preserve = true) {
    02806             // FIXME preserve unimplemented
    02807             BOOST_UBLAS_CHECK (!preserve, internal_logic ());
    02808             size1_ = size1;
    02809             size2_ = size2;
    02810             capacity_ = restrict_capacity (capacity_);
    02811             filled1_ = 1;
    02812             filled2_ = 0;
    02813             index1_data_.resize (layout_type::size_M (size1_, size2_) + 1);
    02814             index2_data_.resize (capacity_);
    02815             value_data_.resize (capacity_);
    02816             index1_data_ [filled1_ - 1] = k_based (filled2_);
    02817             storage_invariants ();
    02818         }
    02819
    02820         // Reserving
    02821         BOOST_UBLAS_INLINE
    02822         void reserve (size_type non_zeros, bool preserve = true) {
    02823             capacity_ = restrict_capacity (non_zeros);
    02824             if (preserve) {
    02825                 index2_data_.resize (capacity_, size_type ());
    02826                 value_data_.resize (capacity_, value_type ());
    02827                 filled2_ = (std::min) (capacity_, filled2_);
    02828             }
    02829             else {
    02830                 index2_data_.resize (capacity_);
    02831                 value_data_.resize (capacity_);
    02832                 filled1_ = 1;
    02833                 filled2_ = 0;
    02834                 index1_data_ [filled1_ - 1] = k_based (filled2_);
    02835             }
    02836             storage_invariants ();
    02837        }
    02838
    02839         // Element support
    02840         BOOST_UBLAS_INLINE
    02841         pointer find_element (size_type i, size_type j) {
    02842             return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i, j));
    02843         }
    02844         BOOST_UBLAS_INLINE
    02845         const_pointer find_element (size_type i, size_type j) const {
    02846             size_type element1 (layout_type::index_M (i, j));
    02847             size_type element2 (layout_type::index_m (i, j));
    02848             if (filled1_ <= element1 + 1)
    02849                 return 0;
    02850             vector_const_subiterator_type itv (index1_data_.begin () + element1);
    02851             const_subiterator_type it_begin (index2_data_.begin () + zero_based (*itv));
    02852             const_subiterator_type it_end (index2_data_.begin () + zero_based (*(itv + 1)));
    02853             const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (element2), std::less<size_type> ()));
    02854             if (it == it_end || *it != k_based (element2))
    02855                 return 0;
    02856             return &value_data_ [it - index2_data_.begin ()];
    02857         }
    02858
    02859         // Element access
    02860         BOOST_UBLAS_INLINE
    02861         const_reference operator () (size_type i, size_type j) const {
    02862             const_pointer p = find_element (i, j);
    02863             if (p)
    02864                 return *p;
    02865             else
    02866                 return zero_;
    02867         }
    02868         BOOST_UBLAS_INLINE
    02869         reference operator () (size_type i, size_type j) {
    02870 #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
    02871             size_type element1 (layout_type::index_M (i, j));
    02872             size_type element2 (layout_type::index_m (i, j));
    02873             if (filled1_ <= element1 + 1)
    02874                 return insert_element (i, j, value_type/*zero*/());
    02875             pointer p = find_element (i, j);
    02876             if (p)
    02877                 return *p;
    02878             else
    02879                 return insert_element (i, j, value_type/*zero*/());
    02880 #else
    02881             return reference (*this, i, j);
    02882 #endif
    02883         }
    02884
    02885         // Element assignment
    02886         BOOST_UBLAS_INLINE
    02887         true_reference insert_element (size_type i, size_type j, const_reference t) {
    02888             BOOST_UBLAS_CHECK (!find_element (i, j), bad_index ());        // duplicate element
    02889             if (filled2_ >= capacity_)
    02890                 reserve (2 * filled2_, true);
    02891             BOOST_UBLAS_CHECK (filled2_ < capacity_, internal_logic ());
    02892             size_type element1 = layout_type::index_M (i, j);
    02893             size_type element2 = layout_type::index_m (i, j);
    02894             while (filled1_ <= element1 + 1) {
    02895                 index1_data_ [filled1_] = k_based (filled2_);
    02896                 ++ filled1_;
    02897             }
    02898             vector_subiterator_type itv (index1_data_.begin () + element1);
    02899             subiterator_type it_begin (index2_data_.begin () + zero_based (*itv));
    02900             subiterator_type it_end (index2_data_.begin () + zero_based (*(itv + 1)));
    02901             subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (element2), std::less<size_type> ()));
    02902             typename std::iterator_traits<subiterator_type>::difference_type n = it - index2_data_.begin ();
    02903             BOOST_UBLAS_CHECK (it == it_end || *it != k_based (element2), internal_logic ());   // duplicate bound by lower_bound
    02904             ++ filled2_;
    02905             it = index2_data_.begin () + n;
    02906             std::copy_backward (it, index2_data_.begin () + filled2_ - 1, index2_data_.begin () + filled2_);
    02907             *it = k_based (element2);
    02908             typename value_array_type::iterator itt (value_data_.begin () + n);
    02909             std::copy_backward (itt, value_data_.begin () + filled2_ - 1, value_data_.begin () + filled2_);
    02910             *itt = t;
    02911             while (element1 + 1 < filled1_) {
    02912                 ++ index1_data_ [element1 + 1];
    02913                 ++ element1;
    02914             }
    02915             storage_invariants ();
    02916             return *itt;
    02917         }
    02918         BOOST_UBLAS_INLINE
    02919         void erase_element (size_type i, size_type j) {
    02920             size_type element1 = layout_type::index_M (i, j);
    02921             size_type element2 = layout_type::index_m (i, j);
    02922             if (element1 + 1 >= filled1_)
    02923                 return;
    02924             vector_subiterator_type itv (index1_data_.begin () + element1);
    02925             subiterator_type it_begin (index2_data_.begin () + zero_based (*itv));
    02926             subiterator_type it_end (index2_data_.begin () + zero_based (*(itv + 1)));
    02927             subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (element2), std::less<size_type> ()));
    02928             if (it != it_end && *it == k_based (element2)) {
    02929                 typename std::iterator_traits<subiterator_type>::difference_type n = it - index2_data_.begin ();
    02930                 std::copy (it + 1, index2_data_.begin () + filled2_, it);
    02931                 typename value_array_type::iterator itt (value_data_.begin () + n);
    02932                 std::copy (itt + 1, value_data_.begin () + filled2_, itt);
    02933                 -- filled2_;
    02934                 while (index1_data_ [filled1_ - 2] > k_based (filled2_)) {
    02935                     index1_data_ [filled1_ - 1] = 0;
    02936                     -- filled1_;
    02937                 }
    02938                 while (element1 + 1 < filled1_) {
    02939                     -- index1_data_ [element1 + 1];
    02940                     ++ element1;
    02941                 }
    02942             }
    02943             storage_invariants ();
    02944         }
    02945
    02946         // Zeroing
    02947         BOOST_UBLAS_INLINE
    02948         void clear () {
    02949             filled1_ = 1;
    02950             filled2_ = 0;
    02951             index1_data_ [filled1_ - 1] = k_based (filled2_);
    02952             storage_invariants ();
    02953         }
    02954
    02955         // Assignment
    02956         BOOST_UBLAS_INLINE
    02957         compressed_matrix &operator = (const compressed_matrix &m) {
    02958             if (this != &m) {
    02959                 size1_ = m.size1_;
    02960                 size2_ = m.size2_;
    02961                 capacity_ = m.capacity_;
    02962                 filled1_ = m.filled1_;
    02963                 filled2_ = m.filled2_;
    02964                 index1_data_ = m.index1_data_;
    02965                 index2_data_ = m.index2_data_;
    02966                 value_data_ = m.value_data_;
    02967             }
    02968             storage_invariants ();
    02969             return *this;
    02970         }
    02971         template<class C>          // Container assignment without temporary
    02972         BOOST_UBLAS_INLINE
    02973         compressed_matrix &operator = (const matrix_container<C> &m) {
    02974             resize (m ().size1 (), m ().size2 (), false);
    02975             assign (m);
    02976             return *this;
    02977         }
    02978         BOOST_UBLAS_INLINE
    02979         compressed_matrix &assign_temporary (compressed_matrix &m) {
    02980             swap (m);
    02981             return *this;
    02982         }
    02983         template<class AE>
    02984         BOOST_UBLAS_INLINE
    02985         compressed_matrix &operator = (const matrix_expression<AE> &ae) {
    02986             self_type temporary (ae, capacity_);
    02987             return assign_temporary (temporary);
    02988         }
    02989         template<class AE>
    02990         BOOST_UBLAS_INLINE
    02991         compressed_matrix &assign (const matrix_expression<AE> &ae) {
    02992             matrix_assign<scalar_assign> (*this, ae);
    02993             return *this;
    02994         }
    02995         template<class AE>
    02996         BOOST_UBLAS_INLINE
    02997         compressed_matrix& operator += (const matrix_expression<AE> &ae) {
    02998             self_type temporary (*this + ae, capacity_);
    02999             return assign_temporary (temporary);
    03000         }
    03001         template<class C>          // Container assignment without temporary
    03002         BOOST_UBLAS_INLINE
    03003         compressed_matrix &operator += (const matrix_container<C> &m) {
    03004             plus_assign (m);
    03005             return *this;
    03006         }
    03007         template<class AE>
    03008         BOOST_UBLAS_INLINE
    03009         compressed_matrix &plus_assign (const matrix_expression<AE> &ae) {
    03010             matrix_assign<scalar_plus_assign> (*this, ae);
    03011             return *this;
    03012         }
    03013         template<class AE>
    03014         BOOST_UBLAS_INLINE
    03015         compressed_matrix& operator -= (const matrix_expression<AE> &ae) {
    03016             self_type temporary (*this - ae, capacity_);
    03017             return assign_temporary (temporary);
    03018         }
    03019         template<class C>          // Container assignment without temporary
    03020         BOOST_UBLAS_INLINE
    03021         compressed_matrix &operator -= (const matrix_container<C> &m) {
    03022             minus_assign (m);
    03023             return *this;
    03024         }
    03025         template<class AE>
    03026         BOOST_UBLAS_INLINE
    03027         compressed_matrix &minus_assign (const matrix_expression<AE> &ae) {
    03028             matrix_assign<scalar_minus_assign> (*this, ae);
    03029             return *this;
    03030         }
    03031         template<class AT>
    03032         BOOST_UBLAS_INLINE
    03033         compressed_matrix& operator *= (const AT &at) {
    03034             matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
    03035             return *this;
    03036         }
    03037         template<class AT>
    03038         BOOST_UBLAS_INLINE
    03039         compressed_matrix& operator /= (const AT &at) {
    03040             matrix_assign_scalar<scalar_divides_assign> (*this, at);
    03041             return *this;
    03042         }
    03043
    03044         // Swapping
    03045         BOOST_UBLAS_INLINE
    03046         void swap (compressed_matrix &m) {
    03047             if (this != &m) {
    03048                 std::swap (size1_, m.size1_);
    03049                 std::swap (size2_, m.size2_);
    03050                 std::swap (capacity_, m.capacity_);
    03051                 std::swap (filled1_, m.filled1_);
    03052                 std::swap (filled2_, m.filled2_);
    03053                 index1_data_.swap (m.index1_data_);
    03054                 index2_data_.swap (m.index2_data_);
    03055                 value_data_.swap (m.value_data_);
    03056             }
    03057             storage_invariants ();
    03058         }
    03059         BOOST_UBLAS_INLINE
    03060         friend void swap (compressed_matrix &m1, compressed_matrix &m2) {
    03061             m1.swap (m2);
    03062         }
    03063
    03064         // Back element insertion and erasure
    03065         BOOST_UBLAS_INLINE
    03066         void push_back (size_type i, size_type j, const_reference t) {
    03067             if (filled2_ >= capacity_)
    03068                 reserve (2 * filled2_, true);
    03069             BOOST_UBLAS_CHECK (filled2_ < capacity_, internal_logic ());
    03070             size_type element1 = layout_type::index_M (i, j);
    03071             size_type element2 = layout_type::index_m (i, j);
    03072             while (filled1_ < element1 + 2) {
    03073                 index1_data_ [filled1_] = k_based (filled2_);
    03074                 ++ filled1_;
    03075             }
    03076             // must maintain sort order
    03077             BOOST_UBLAS_CHECK ((filled1_ == element1 + 2 &&
    03078                                 (filled2_ == zero_based (index1_data_ [filled1_ - 2]) ||
    03079                                 index2_data_ [filled2_ - 1] < k_based (element2))), external_logic ());
    03080             ++ filled2_;
    03081             index1_data_ [filled1_ - 1] = k_based (filled2_);
    03082             index2_data_ [filled2_ - 1] = k_based (element2);
    03083             value_data_ [filled2_ - 1] = t;
    03084             storage_invariants ();
    03085         }
    03086         BOOST_UBLAS_INLINE
    03087         void pop_back () {
    03088             BOOST_UBLAS_CHECK (filled1_ > 0 && filled2_ > 0, external_logic ());
    03089             -- filled2_;
    03090             while (index1_data_ [filled1_ - 2] > k_based (filled2_)) {
    03091                 index1_data_ [filled1_ - 1] = 0;
    03092                 -- filled1_;
    03093             }
    03094             -- index1_data_ [filled1_ - 1];
    03095             storage_invariants ();
    03096         }
    03097
    03098         // Iterator types
    03099     private:
    03100         // Use index array iterator
    03101         typedef typename IA::const_iterator vector_const_subiterator_type;
    03102         typedef typename IA::iterator vector_subiterator_type;
    03103         typedef typename IA::const_iterator const_subiterator_type;
    03104         typedef typename IA::iterator subiterator_type;
    03105
    03106         BOOST_UBLAS_INLINE
    03107         true_reference at_element (size_type i, size_type j) {
    03108             pointer p = find_element (i, j);
    03109             BOOST_UBLAS_CHECK (p, bad_index ());
    03110             return *p;
    03111         }
    03112
    03113     public:
    03114         class const_iterator1;
    03115         class iterator1;
    03116         class const_iterator2;
    03117         class iterator2;
    03118         typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
    03119         typedef reverse_iterator_base1<iterator1> reverse_iterator1;
    03120         typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
    03121         typedef reverse_iterator_base2<iterator2> reverse_iterator2;
    03122
    03123         // Element lookup
    03124         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    03125         const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const {
    03126             for (;;) {
    03127                 array_size_type address1 (layout_type::index_M (i, j));
    03128                 array_size_type address2 (layout_type::index_m (i, j));
    03129                 vector_const_subiterator_type itv (index1_data_.begin () + (std::min) (filled1_ - 1, address1));
    03130                 if (filled1_ <= address1 + 1)
    03131                     return const_iterator1 (*this, rank, i, j, itv, index2_data_.begin () + filled2_);
    03132
    03133                 const_subiterator_type it_begin (index2_data_.begin () + zero_based (*itv));
    03134                 const_subiterator_type it_end (index2_data_.begin () + zero_based (*(itv + 1)));
    03135
    03136                 const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ()));
    03137                 if (rank == 0)
    03138                     return const_iterator1 (*this, rank, i, j, itv, it);
    03139                 if (it != it_end && zero_based (*it) == address2)
    03140                     return const_iterator1 (*this, rank, i, j, itv, it);
    03141                 if (direction > 0) {
    03142                     if (layout_type::fast_i ()) {
    03143                         if (it == it_end)
    03144                             return const_iterator1 (*this, rank, i, j, itv, it);
    03145                         i = zero_based (*it);
    03146                     } else {
    03147                         if (i >= size1_)
    03148                             return const_iterator1 (*this, rank, i, j, itv, it);
    03149                         ++ i;
    03150                     }
    03151                 } else /* if (direction < 0)  */ {
    03152                     if (layout_type::fast_i ()) {
    03153                         if (it == index2_data_.begin () + zero_based (*itv))
    03154                             return const_iterator1 (*this, rank, i, j, itv, it);
    03155                         i = zero_based (*(it - 1));
    03156                     } else {
    03157                         if (i == 0)
    03158                             return const_iterator1 (*this, rank, i, j, itv, it);
    03159                         -- i;
    03160                     }
    03161                 }
    03162             }
    03163         }
    03164         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    03165         iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) {
    03166             for (;;) {
    03167                 array_size_type address1 (layout_type::index_M (i, j));
    03168                 array_size_type address2 (layout_type::index_m (i, j));
    03169                 vector_subiterator_type itv (index1_data_.begin () + (std::min) (filled1_ - 1, address1));
    03170                 if (filled1_ <= address1 + 1)
    03171                     return iterator1 (*this, rank, i, j, itv, index2_data_.begin () + filled2_);
    03172
    03173                 subiterator_type it_begin (index2_data_.begin () + zero_based (*itv));
    03174                 subiterator_type it_end (index2_data_.begin () + zero_based (*(itv + 1)));
    03175
    03176                 subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ()));
    03177                 if (rank == 0)
    03178                     return iterator1 (*this, rank, i, j, itv, it);
    03179                 if (it != it_end && zero_based (*it) == address2)
    03180                     return iterator1 (*this, rank, i, j, itv, it);
    03181                 if (direction > 0) {
    03182                     if (layout_type::fast_i ()) {
    03183                         if (it == it_end)
    03184                             return iterator1 (*this, rank, i, j, itv, it);
    03185                         i = zero_based (*it);
    03186                     } else {
    03187                         if (i >= size1_)
    03188                             return iterator1 (*this, rank, i, j, itv, it);
    03189                         ++ i;
    03190                     }
    03191                 } else /* if (direction < 0)  */ {
    03192                     if (layout_type::fast_i ()) {
    03193                         if (it == index2_data_.begin () + zero_based (*itv))
    03194                             return iterator1 (*this, rank, i, j, itv, it);
    03195                         i = zero_based (*(it - 1));
    03196                     } else {
    03197                         if (i == 0)
    03198                             return iterator1 (*this, rank, i, j, itv, it);
    03199                         -- i;
    03200                     }
    03201                 }
    03202             }
    03203         }
    03204         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    03205         const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const {
    03206             for (;;) {
    03207                 array_size_type address1 (layout_type::index_M (i, j));
    03208                 array_size_type address2 (layout_type::index_m (i, j));
    03209                 vector_const_subiterator_type itv (index1_data_.begin () + (std::min) (filled1_ - 1, address1));
    03210                 if (filled1_ <= address1 + 1)
    03211                     return const_iterator2 (*this, rank, i, j, itv, index2_data_.begin () + filled2_);
    03212
    03213                 const_subiterator_type it_begin (index2_data_.begin () + zero_based (*itv));
    03214                 const_subiterator_type it_end (index2_data_.begin () + zero_based (*(itv + 1)));
    03215
    03216                 const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ()));
    03217                 if (rank == 0)
    03218                     return const_iterator2 (*this, rank, i, j, itv, it);
    03219                 if (it != it_end && zero_based (*it) == address2)
    03220                     return const_iterator2 (*this, rank, i, j, itv, it);
    03221                 if (direction > 0) {
    03222                     if (layout_type::fast_j ()) {
    03223                         if (it == it_end)
    03224                             return const_iterator2 (*this, rank, i, j, itv, it);
    03225                         j = zero_based (*it);
    03226                     } else {
    03227                         if (j >= size2_)
    03228                             return const_iterator2 (*this, rank, i, j, itv, it);
    03229                         ++ j;
    03230                     }
    03231                 } else /* if (direction < 0)  */ {
    03232                     if (layout_type::fast_j ()) {
    03233                         if (it == index2_data_.begin () + zero_based (*itv))
    03234                             return const_iterator2 (*this, rank, i, j, itv, it);
    03235                         j = zero_based (*(it - 1));
    03236                     } else {
    03237                         if (j == 0)
    03238                             return const_iterator2 (*this, rank, i, j, itv, it);
    03239                         -- j;
    03240                     }
    03241                 }
    03242             }
    03243         }
    03244         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    03245         iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) {
    03246             for (;;) {
    03247                 array_size_type address1 (layout_type::index_M (i, j));
    03248                 array_size_type address2 (layout_type::index_m (i, j));
    03249                 vector_subiterator_type itv (index1_data_.begin () + (std::min) (filled1_ - 1, address1));
    03250                 if (filled1_ <= address1 + 1)
    03251                     return iterator2 (*this, rank, i, j, itv, index2_data_.begin () + filled2_);
    03252
    03253                 subiterator_type it_begin (index2_data_.begin () + zero_based (*itv));
    03254                 subiterator_type it_end (index2_data_.begin () + zero_based (*(itv + 1)));
    03255
    03256                 subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ()));
    03257                 if (rank == 0)
    03258                     return iterator2 (*this, rank, i, j, itv, it);
    03259                 if (it != it_end && zero_based (*it) == address2)
    03260                     return iterator2 (*this, rank, i, j, itv, it);
    03261                 if (direction > 0) {
    03262                     if (layout_type::fast_j ()) {
    03263                         if (it == it_end)
    03264                             return iterator2 (*this, rank, i, j, itv, it);
    03265                         j = zero_based (*it);
    03266                     } else {
    03267                         if (j >= size2_)
    03268                             return iterator2 (*this, rank, i, j, itv, it);
    03269                         ++ j;
    03270                     }
    03271                 } else /* if (direction < 0)  */ {
    03272                     if (layout_type::fast_j ()) {
    03273                         if (it == index2_data_.begin () + zero_based (*itv))
    03274                             return iterator2 (*this, rank, i, j, itv, it);
    03275                         j = zero_based (*(it - 1));
    03276                     } else {
    03277                         if (j == 0)
    03278                             return iterator2 (*this, rank, i, j, itv, it);
    03279                         -- j;
    03280                     }
    03281                 }
    03282             }
    03283         }
    03284
    03285
    03286         class const_iterator1:
    03287             public container_const_reference<compressed_matrix>,
    03288             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
    03289                                                const_iterator1, value_type> {
    03290         public:
    03291             typedef typename compressed_matrix::value_type value_type;
    03292             typedef typename compressed_matrix::difference_type difference_type;
    03293             typedef typename compressed_matrix::const_reference reference;
    03294             typedef const typename compressed_matrix::pointer pointer;
    03295
    03296             typedef const_iterator2 dual_iterator_type;
    03297             typedef const_reverse_iterator2 dual_reverse_iterator_type;
    03298
    03299             // Construction and destruction
    03300             BOOST_UBLAS_INLINE
    03301             const_iterator1 ():
    03302                 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
    03303             BOOST_UBLAS_INLINE
    03304             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):
    03305                 container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
    03306             BOOST_UBLAS_INLINE
    03307             const_iterator1 (const iterator1 &it):
    03308                 container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {}
    03309
    03310             // Arithmetic
    03311             BOOST_UBLAS_INLINE
    03312             const_iterator1 &operator ++ () {
    03313                 if (rank_ == 1 && layout_type::fast_i ())
    03314                     ++ it_;
    03315                 else {
    03316                     i_ = index1 () + 1;
    03317                     if (rank_ == 1)
    03318                         *this = (*this) ().find1 (rank_, i_, j_, 1);
    03319                 }
    03320                 return *this;
    03321             }
    03322             BOOST_UBLAS_INLINE
    03323             const_iterator1 &operator -- () {
    03324                 if (rank_ == 1 && layout_type::fast_i ())
    03325                     -- it_;
    03326                 else {
    03327                     --i_;
    03328                     if (rank_ == 1)
    03329                         *this = (*this) ().find1 (rank_, i_, j_, -1);
    03330                 }
    03331                 return *this;
    03332             }
    03333
    03334             // Dereference
    03335             BOOST_UBLAS_INLINE
    03336             const_reference operator * () const {
    03337                 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
    03338                 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
    03339                 if (rank_ == 1) {
    03340                     return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
    03341                 } else {
    03342                     return (*this) () (i_, j_);
    03343                 }
    03344             }
    03345
    03346 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    03347             BOOST_UBLAS_INLINE
    03348 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    03349             typename self_type::
    03350 #endif
    03351             const_iterator2 begin () const {
    03352                 const self_type &m = (*this) ();
    03353                 return m.find2 (1, index1 (), 0);
    03354             }
    03355             BOOST_UBLAS_INLINE
    03356 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    03357             typename self_type::
    03358 #endif
    03359             const_iterator2 end () const {
    03360                 const self_type &m = (*this) ();
    03361                 return m.find2 (1, index1 (), m.size2 ());
    03362             }
    03363             BOOST_UBLAS_INLINE
    03364 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    03365             typename self_type::
    03366 #endif
    03367             const_reverse_iterator2 rbegin () const {
    03368                 return const_reverse_iterator2 (end ());
    03369             }
    03370             BOOST_UBLAS_INLINE
    03371 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    03372             typename self_type::
    03373 #endif
    03374             const_reverse_iterator2 rend () const {
    03375                 return const_reverse_iterator2 (begin ());
    03376             }
    03377 #endif
    03378
    03379             // Indices
    03380             BOOST_UBLAS_INLINE
    03381             size_type index1 () const {
    03382                 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
    03383                 if (rank_ == 1) {
    03384                     BOOST_UBLAS_CHECK (layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
    03385                     return layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
    03386                 } else {
    03387                     return i_;
    03388                 }
    03389             }
    03390             BOOST_UBLAS_INLINE
    03391             size_type index2 () const {
    03392                 if (rank_ == 1) {
    03393                     BOOST_UBLAS_CHECK (layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
    03394                     return layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
    03395                 } else {
    03396                     return j_;
    03397                 }
    03398             }
    03399
    03400             // Assignment
    03401             BOOST_UBLAS_INLINE
    03402             const_iterator1 &operator = (const const_iterator1 &it) {
    03403                 container_const_reference<self_type>::assign (&it ());
    03404                 rank_ = it.rank_;
    03405                 i_ = it.i_;
    03406                 j_ = it.j_;
    03407                 itv_ = it.itv_;
    03408                 it_ = it.it_;
    03409                 return *this;
    03410             }
    03411
    03412             // Comparison
    03413             BOOST_UBLAS_INLINE
    03414             bool operator == (const const_iterator1 &it) const {
    03415                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
    03416                 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
    03417                 if (rank_ == 1 || it.rank_ == 1) {
    03418                     return it_ == it.it_;
    03419                 } else {
    03420                     return i_ == it.i_ && j_ == it.j_;
    03421                 }
    03422             }
    03423
    03424         private:
    03425             int rank_;
    03426             size_type i_;
    03427             size_type j_;
    03428             vector_const_subiterator_type itv_;
    03429             const_subiterator_type it_;
    03430         };
    03431
    03432         BOOST_UBLAS_INLINE
    03433         const_iterator1 begin1 () const {
    03434             return find1 (0, 0, 0);
    03435         }
    03436         BOOST_UBLAS_INLINE
    03437         const_iterator1 end1 () const {
    03438             return find1 (0, size1_, 0);
    03439         }
    03440
    03441         class iterator1:
    03442             public container_reference<compressed_matrix>,
    03443             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
    03444                                                iterator1, value_type> {
    03445         public:
    03446             typedef typename compressed_matrix::value_type value_type;
    03447             typedef typename compressed_matrix::difference_type difference_type;
    03448             typedef typename compressed_matrix::true_reference reference;
    03449             typedef typename compressed_matrix::pointer pointer;
    03450
    03451             typedef iterator2 dual_iterator_type;
    03452             typedef reverse_iterator2 dual_reverse_iterator_type;
    03453
    03454             // Construction and destruction
    03455             BOOST_UBLAS_INLINE
    03456             iterator1 ():
    03457                 container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
    03458             BOOST_UBLAS_INLINE
    03459             iterator1 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it):
    03460                 container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
    03461
    03462             // Arithmetic
    03463             BOOST_UBLAS_INLINE
    03464             iterator1 &operator ++ () {
    03465                 if (rank_ == 1 && layout_type::fast_i ())
    03466                     ++ it_;
    03467                 else {
    03468                     i_ = index1 () + 1;
    03469                     if (rank_ == 1)
    03470                         *this = (*this) ().find1 (rank_, i_, j_, 1);
    03471                 }
    03472                 return *this;
    03473             }
    03474             BOOST_UBLAS_INLINE
    03475             iterator1 &operator -- () {
    03476                 if (rank_ == 1 && layout_type::fast_i ())
    03477                     -- it_;
    03478                 else {
    03479                     --i_;
    03480                     if (rank_ == 1)
    03481                         *this = (*this) ().find1 (rank_, i_, j_, -1);
    03482                 }
    03483                 return *this;
    03484             }
    03485
    03486             // Dereference
    03487             BOOST_UBLAS_INLINE
    03488             reference operator * () const {
    03489                 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
    03490                 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
    03491                 if (rank_ == 1) {
    03492                     return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
    03493                 } else {
    03494                     return (*this) ().at_element (i_, j_);
    03495                 }
    03496             }
    03497
    03498 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    03499             BOOST_UBLAS_INLINE
    03500 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    03501             typename self_type::
    03502 #endif
    03503             iterator2 begin () const {
    03504                 self_type &m = (*this) ();
    03505                 return m.find2 (1, index1 (), 0);
    03506             }
    03507             BOOST_UBLAS_INLINE
    03508 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    03509             typename self_type::
    03510 #endif
    03511             iterator2 end () const {
    03512                 self_type &m = (*this) ();
    03513                 return m.find2 (1, index1 (), m.size2 ());
    03514             }
    03515             BOOST_UBLAS_INLINE
    03516 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    03517             typename self_type::
    03518 #endif
    03519             reverse_iterator2 rbegin () const {
    03520                 return reverse_iterator2 (end ());
    03521             }
    03522             BOOST_UBLAS_INLINE
    03523 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    03524             typename self_type::
    03525 #endif
    03526             reverse_iterator2 rend () const {
    03527                 return reverse_iterator2 (begin ());
    03528             }
    03529 #endif
    03530
    03531             // Indices
    03532             BOOST_UBLAS_INLINE
    03533             size_type index1 () const {
    03534                 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
    03535                 if (rank_ == 1) {
    03536                     BOOST_UBLAS_CHECK (layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
    03537                     return layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
    03538                 } else {
    03539                     return i_;
    03540                 }
    03541             }
    03542             BOOST_UBLAS_INLINE
    03543             size_type index2 () const {
    03544                 if (rank_ == 1) {
    03545                     BOOST_UBLAS_CHECK (layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
    03546                     return layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
    03547                 } else {
    03548                     return j_;
    03549                 }
    03550             }
    03551
    03552             // Assignment
    03553             BOOST_UBLAS_INLINE
    03554             iterator1 &operator = (const iterator1 &it) {
    03555                 container_reference<self_type>::assign (&it ());
    03556                 rank_ = it.rank_;
    03557                 i_ = it.i_;
    03558                 j_ = it.j_;
    03559                 itv_ = it.itv_;
    03560                 it_ = it.it_;
    03561                 return *this;
    03562             }
    03563
    03564             // Comparison
    03565             BOOST_UBLAS_INLINE
    03566             bool operator == (const iterator1 &it) const {
    03567                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
    03568                 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
    03569                 if (rank_ == 1 || it.rank_ == 1) {
    03570                     return it_ == it.it_;
    03571                 } else {
    03572                     return i_ == it.i_ && j_ == it.j_;
    03573                 }
    03574             }
    03575
    03576         private:
    03577             int rank_;
    03578             size_type i_;
    03579             size_type j_;
    03580             vector_subiterator_type itv_;
    03581             subiterator_type it_;
    03582
    03583             friend class const_iterator1;
    03584         };
    03585
    03586         BOOST_UBLAS_INLINE
    03587         iterator1 begin1 () {
    03588             return find1 (0, 0, 0);
    03589         }
    03590         BOOST_UBLAS_INLINE
    03591         iterator1 end1 () {
    03592             return find1 (0, size1_, 0);
    03593         }
    03594
    03595         class const_iterator2:
    03596             public container_const_reference<compressed_matrix>,
    03597             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
    03598                                                const_iterator2, value_type> {
    03599         public:
    03600             typedef typename compressed_matrix::value_type value_type;
    03601             typedef typename compressed_matrix::difference_type difference_type;
    03602             typedef typename compressed_matrix::const_reference reference;
    03603             typedef const typename compressed_matrix::pointer pointer;
    03604
    03605             typedef const_iterator1 dual_iterator_type;
    03606             typedef const_reverse_iterator1 dual_reverse_iterator_type;
    03607
    03608             // Construction and destruction
    03609             BOOST_UBLAS_INLINE
    03610             const_iterator2 ():
    03611                 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
    03612             BOOST_UBLAS_INLINE
    03613             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):
    03614                 container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
    03615             BOOST_UBLAS_INLINE
    03616             const_iterator2 (const iterator2 &it):
    03617                 container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {}
    03618
    03619             // Arithmetic
    03620             BOOST_UBLAS_INLINE
    03621             const_iterator2 &operator ++ () {
    03622                 if (rank_ == 1 && layout_type::fast_j ())
    03623                     ++ it_;
    03624                 else {
    03625                     j_ = index2 () + 1;
    03626                     if (rank_ == 1)
    03627                         *this = (*this) ().find2 (rank_, i_, j_, 1);
    03628                 }
    03629                 return *this;
    03630             }
    03631             BOOST_UBLAS_INLINE
    03632             const_iterator2 &operator -- () {
    03633                 if (rank_ == 1 && layout_type::fast_j ())
    03634                     -- it_;
    03635                 else {
    03636                     --j_;
    03637                     if (rank_ == 1)
    03638                         *this = (*this) ().find2 (rank_, i_, j_, -1);
    03639                 }
    03640                 return *this;
    03641             }
    03642
    03643             // Dereference
    03644             BOOST_UBLAS_INLINE
    03645             const_reference operator * () const {
    03646                 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
    03647                 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
    03648                 if (rank_ == 1) {
    03649                     return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
    03650                 } else {
    03651                     return (*this) () (i_, j_);
    03652                 }
    03653             }
    03654
    03655 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    03656             BOOST_UBLAS_INLINE
    03657 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    03658             typename self_type::
    03659 #endif
    03660             const_iterator1 begin () const {
    03661                 const self_type &m = (*this) ();
    03662                 return m.find1 (1, 0, index2 ());
    03663             }
    03664             BOOST_UBLAS_INLINE
    03665 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    03666             typename self_type::
    03667 #endif
    03668             const_iterator1 end () const {
    03669                 const self_type &m = (*this) ();
    03670                 return m.find1 (1, m.size1 (), index2 ());
    03671             }
    03672             BOOST_UBLAS_INLINE
    03673 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    03674             typename self_type::
    03675 #endif
    03676             const_reverse_iterator1 rbegin () const {
    03677                 return const_reverse_iterator1 (end ());
    03678             }
    03679             BOOST_UBLAS_INLINE
    03680 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    03681             typename self_type::
    03682 #endif
    03683             const_reverse_iterator1 rend () const {
    03684                 return const_reverse_iterator1 (begin ());
    03685             }
    03686 #endif
    03687
    03688             // Indices
    03689             BOOST_UBLAS_INLINE
    03690             size_type index1 () const {
    03691                 if (rank_ == 1) {
    03692                     BOOST_UBLAS_CHECK (layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
    03693                     return layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
    03694                 } else {
    03695                     return i_;
    03696                 }
    03697             }
    03698             BOOST_UBLAS_INLINE
    03699             size_type index2 () const {
    03700                 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
    03701                 if (rank_ == 1) {
    03702                     BOOST_UBLAS_CHECK (layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
    03703                     return layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
    03704                 } else {
    03705                     return j_;
    03706                 }
    03707             }
    03708
    03709             // Assignment
    03710             BOOST_UBLAS_INLINE
    03711             const_iterator2 &operator = (const const_iterator2 &it) {
    03712                 container_const_reference<self_type>::assign (&it ());
    03713                 rank_ = it.rank_;
    03714                 i_ = it.i_;
    03715                 j_ = it.j_;
    03716                 itv_ = it.itv_;
    03717                 it_ = it.it_;
    03718                 return *this;
    03719             }
    03720
    03721             // Comparison
    03722             BOOST_UBLAS_INLINE
    03723             bool operator == (const const_iterator2 &it) const {
    03724                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
    03725                 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
    03726                 if (rank_ == 1 || it.rank_ == 1) {
    03727                     return it_ == it.it_;
    03728                 } else {
    03729                     return i_ == it.i_ && j_ == it.j_;
    03730                 }
    03731             }
    03732
    03733         private:
    03734             int rank_;
    03735             size_type i_;
    03736             size_type j_;
    03737             vector_const_subiterator_type itv_;
    03738             const_subiterator_type it_;
    03739         };
    03740
    03741         BOOST_UBLAS_INLINE
    03742         const_iterator2 begin2 () const {
    03743             return find2 (0, 0, 0);
    03744         }
    03745         BOOST_UBLAS_INLINE
    03746         const_iterator2 end2 () const {
    03747             return find2 (0, 0, size2_);
    03748         }
    03749
    03750         class iterator2:
    03751             public container_reference<compressed_matrix>,
    03752             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
    03753                                                iterator2, value_type> {
    03754         public:
    03755             typedef typename compressed_matrix::value_type value_type;
    03756             typedef typename compressed_matrix::difference_type difference_type;
    03757             typedef typename compressed_matrix::true_reference reference;
    03758             typedef typename compressed_matrix::pointer pointer;
    03759
    03760             typedef iterator1 dual_iterator_type;
    03761             typedef reverse_iterator1 dual_reverse_iterator_type;
    03762
    03763             // Construction and destruction
    03764             BOOST_UBLAS_INLINE
    03765             iterator2 ():
    03766                 container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
    03767             BOOST_UBLAS_INLINE
    03768             iterator2 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it):
    03769                 container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
    03770
    03771             // Arithmetic
    03772             BOOST_UBLAS_INLINE
    03773             iterator2 &operator ++ () {
    03774                 if (rank_ == 1 && layout_type::fast_j ())
    03775                     ++ it_;
    03776                 else {
    03777                     j_ = index2 () + 1;
    03778                     if (rank_ == 1)
    03779                         *this = (*this) ().find2 (rank_, i_, j_, 1);
    03780                 }
    03781                 return *this;
    03782             }
    03783             BOOST_UBLAS_INLINE
    03784             iterator2 &operator -- () {
    03785                 if (rank_ == 1 && layout_type::fast_j ())
    03786                     -- it_;
    03787                 else {
    03788                     --j_;
    03789                     if (rank_ == 1)
    03790                         *this = (*this) ().find2 (rank_, i_, j_, -1);
    03791                 }
    03792                 return *this;
    03793             }
    03794
    03795             // Dereference
    03796             BOOST_UBLAS_INLINE
    03797             reference operator * () const {
    03798                 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
    03799                 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
    03800                 if (rank_ == 1) {
    03801                     return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
    03802                 } else {
    03803                     return (*this) ().at_element (i_, j_);
    03804                 }
    03805             }
    03806
    03807 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    03808             BOOST_UBLAS_INLINE
    03809 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    03810             typename self_type::
    03811 #endif
    03812             iterator1 begin () const {
    03813                 self_type &m = (*this) ();
    03814                 return m.find1 (1, 0, index2 ());
    03815             }
    03816             BOOST_UBLAS_INLINE
    03817 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    03818             typename self_type::
    03819 #endif
    03820             iterator1 end () const {
    03821                 self_type &m = (*this) ();
    03822                 return m.find1 (1, m.size1 (), index2 ());
    03823             }
    03824             BOOST_UBLAS_INLINE
    03825 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    03826             typename self_type::
    03827 #endif
    03828             reverse_iterator1 rbegin () const {
    03829                 return reverse_iterator1 (end ());
    03830             }
    03831             BOOST_UBLAS_INLINE
    03832 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    03833             typename self_type::
    03834 #endif
    03835             reverse_iterator1 rend () const {
    03836                 return reverse_iterator1 (begin ());
    03837             }
    03838 #endif
    03839
    03840             // Indices
    03841             BOOST_UBLAS_INLINE
    03842             size_type index1 () const {
    03843                 if (rank_ == 1) {
    03844                     BOOST_UBLAS_CHECK (layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
    03845                     return layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
    03846                 } else {
    03847                     return i_;
    03848                 }
    03849             }
    03850             BOOST_UBLAS_INLINE
    03851             size_type index2 () const {
    03852                 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
    03853                 if (rank_ == 1) {
    03854                     BOOST_UBLAS_CHECK (layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
    03855                     return layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
    03856                 } else {
    03857                     return j_;
    03858                 }
    03859             }
    03860
    03861             // Assignment
    03862             BOOST_UBLAS_INLINE
    03863             iterator2 &operator = (const iterator2 &it) {
    03864                 container_reference<self_type>::assign (&it ());
    03865                 rank_ = it.rank_;
    03866                 i_ = it.i_;
    03867                 j_ = it.j_;
    03868                 itv_ = it.itv_;
    03869                 it_ = it.it_;
    03870                 return *this;
    03871             }
    03872
    03873             // Comparison
    03874             BOOST_UBLAS_INLINE
    03875             bool operator == (const iterator2 &it) const {
    03876                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
    03877                 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
    03878                 if (rank_ == 1 || it.rank_ == 1) {
    03879                     return it_ == it.it_;
    03880                 } else {
    03881                     return i_ == it.i_ && j_ == it.j_;
    03882                 }
    03883             }
    03884
    03885         private:
    03886             int rank_;
    03887             size_type i_;
    03888             size_type j_;
    03889             vector_subiterator_type itv_;
    03890             subiterator_type it_;
    03891
    03892             friend class const_iterator2;
    03893         };
    03894
    03895         BOOST_UBLAS_INLINE
    03896         iterator2 begin2 () {
    03897             return find2 (0, 0, 0);
    03898         }
    03899         BOOST_UBLAS_INLINE
    03900         iterator2 end2 () {
    03901             return find2 (0, 0, size2_);
    03902         }
    03903
    03904         // Reverse iterators
    03905
    03906         BOOST_UBLAS_INLINE
    03907         const_reverse_iterator1 rbegin1 () const {
    03908             return const_reverse_iterator1 (end1 ());
    03909         }
    03910         BOOST_UBLAS_INLINE
    03911         const_reverse_iterator1 rend1 () const {
    03912             return const_reverse_iterator1 (begin1 ());
    03913         }
    03914
    03915         BOOST_UBLAS_INLINE
    03916         reverse_iterator1 rbegin1 () {
    03917             return reverse_iterator1 (end1 ());
    03918         }
    03919         BOOST_UBLAS_INLINE
    03920         reverse_iterator1 rend1 () {
    03921             return reverse_iterator1 (begin1 ());
    03922         }
    03923
    03924         BOOST_UBLAS_INLINE
    03925         const_reverse_iterator2 rbegin2 () const {
    03926             return const_reverse_iterator2 (end2 ());
    03927         }
    03928         BOOST_UBLAS_INLINE
    03929         const_reverse_iterator2 rend2 () const {
    03930             return const_reverse_iterator2 (begin2 ());
    03931         }
    03932
    03933         BOOST_UBLAS_INLINE
    03934         reverse_iterator2 rbegin2 () {
    03935             return reverse_iterator2 (end2 ());
    03936         }
    03937         BOOST_UBLAS_INLINE
    03938         reverse_iterator2 rend2 () {
    03939             return reverse_iterator2 (begin2 ());
    03940         }
    03941
    03942          // Serialization
    03943         template<class Archive>
    03944         void serialize(Archive & ar, const unsigned int /* file_version */){
    03945             serialization::collection_size_type s1 (size1_);
    03946             serialization::collection_size_type s2 (size2_);
    03947             ar & serialization::make_nvp("size1",s1);
    03948             ar & serialization::make_nvp("size2",s2);
    03949             if (Archive::is_loading::value) {
    03950                 size1_ = s1;
    03951                 size2_ = s2;
    03952             }
    03953             ar & serialization::make_nvp("capacity", capacity_);
    03954             ar & serialization::make_nvp("filled1", filled1_);
    03955             ar & serialization::make_nvp("filled2", filled2_);
    03956             ar & serialization::make_nvp("index1_data", index1_data_);
    03957             ar & serialization::make_nvp("index2_data", index2_data_);
    03958             ar & serialization::make_nvp("value_data", value_data_);
    03959             storage_invariants();
    03960         }
    03961
    03962     private:
    03963         void storage_invariants () const {
    03964             BOOST_UBLAS_CHECK (layout_type::size_M (size1_, size2_) + 1 == index1_data_.size (), internal_logic ());
    03965             BOOST_UBLAS_CHECK (capacity_ == index2_data_.size (), internal_logic ());
    03966             BOOST_UBLAS_CHECK (capacity_ == value_data_.size (), internal_logic ());
    03967             BOOST_UBLAS_CHECK (filled1_ > 0 && filled1_ <= layout_type::size_M (size1_, size2_) + 1, internal_logic ());
    03968             BOOST_UBLAS_CHECK (filled2_ <= capacity_, internal_logic ());
    03969             BOOST_UBLAS_CHECK (index1_data_ [filled1_ - 1] == k_based (filled2_), internal_logic ());
    03970         }
    03971
    03972         size_type size1_;
    03973         size_type size2_;
    03974         array_size_type capacity_;
    03975         array_size_type filled1_;
    03976         array_size_type filled2_;
    03977         index_array_type index1_data_;
    03978         index_array_type index2_data_;
    03979         value_array_type value_data_;
    03980         static const value_type zero_;
    03981
    03982         BOOST_UBLAS_INLINE
    03983         static size_type zero_based (size_type k_based_index) {
    03984             return k_based_index - IB;
    03985         }
    03986         BOOST_UBLAS_INLINE
    03987         static size_type k_based (size_type zero_based_index) {
    03988             return zero_based_index + IB;
    03989         }
    03990
    03991         friend class iterator1;
    03992         friend class iterator2;
    03993         friend class const_iterator1;
    03994         friend class const_iterator2;
    03995     };
    03996
    03997     template<class T, class L, std::size_t IB, class IA, class TA>
    03998     const typename compressed_matrix<T, L, IB, IA, TA>::value_type compressed_matrix<T, L, IB, IA, TA>::zero_ = value_type/*zero*/();
    03999
    04000
    04001     // Coordinate array based sparse matrix class
    04002     // Thanks to Kresimir Fresl for extending this to cover different index bases.
    04003     template<class T, class L, std::size_t IB, class IA, class TA>
    04004     class coordinate_matrix:
    04005         public matrix_container<coordinate_matrix<T, L, IB, IA, TA> > {
    04006
    04007         typedef T &true_reference;
    04008         typedef T *pointer;
    04009         typedef const T *const_pointer;
    04010         typedef L layout_type;
    04011         typedef coordinate_matrix<T, L, IB, IA, TA> self_type;
    04012     public:
    04013 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
    04014         using matrix_container<self_type>::operator ();
    04015 #endif
    04016         // ISSUE require type consistency check, is_convertable (IA::size_type, TA::size_type)
    04017         typedef typename IA::value_type size_type;
    04018         // ISSUE difference_type cannot be deduced for sparse indices, we only know the value_type
    04019         typedef std::ptrdiff_t difference_type;
    04020         // size_type for the data arrays.
    04021         typedef typename IA::size_type array_size_type;
    04022         typedef T value_type;
    04023         typedef const T &const_reference;
    04024 #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
    04025         typedef T &reference;
    04026 #else
    04027         typedef sparse_matrix_element<self_type> reference;
    04028 #endif
    04029         typedef IA index_array_type;
    04030         typedef TA value_array_type;
    04031         typedef const matrix_reference<const self_type> const_closure_type;
    04032         typedef matrix_reference<self_type> closure_type;
    04033         typedef coordinate_vector<T, IB, IA, TA> vector_temporary_type;
    04034         typedef self_type matrix_temporary_type;
    04035         typedef sparse_tag storage_category;
    04036         typedef typename L::orientation_category orientation_category;
    04037
    04038         // Construction and destruction
    04039         BOOST_UBLAS_INLINE
    04040         coordinate_matrix ():
    04041             matrix_container<self_type> (),
    04042             size1_ (0), size2_ (0), capacity_ (restrict_capacity (0)),
    04043             filled_ (0), sorted_filled_ (filled_), sorted_ (true),
    04044             index1_data_ (capacity_), index2_data_ (capacity_), value_data_ (capacity_) {
    04045             storage_invariants ();
    04046         }
    04047         BOOST_UBLAS_INLINE
    04048         coordinate_matrix (size_type size1, size_type size2, array_size_type non_zeros = 0):
    04049             matrix_container<self_type> (),
    04050             size1_ (size1), size2_ (size2), capacity_ (restrict_capacity (non_zeros)),
    04051             filled_ (0), sorted_filled_ (filled_), sorted_ (true),
    04052             index1_data_ (capacity_), index2_data_ (capacity_), value_data_ (capacity_) {
    04053             storage_invariants ();
    04054         }
    04055         BOOST_UBLAS_INLINE
    04056         coordinate_matrix (const coordinate_matrix &m):
    04057             matrix_container<self_type> (),
    04058             size1_ (m.size1_), size2_ (m.size2_), capacity_ (m.capacity_),
    04059             filled_ (m.filled_), sorted_filled_ (m.sorted_filled_), sorted_ (m.sorted_),
    04060             index1_data_ (m.index1_data_), index2_data_ (m.index2_data_), value_data_ (m.value_data_) {
    04061             storage_invariants ();
    04062         }
    04063         template<class AE>
    04064         BOOST_UBLAS_INLINE
    04065         coordinate_matrix (const matrix_expression<AE> &ae, array_size_type non_zeros = 0):
    04066             matrix_container<self_type> (),
    04067             size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), capacity_ (restrict_capacity (non_zeros)),
    04068             filled_ (0), sorted_filled_ (filled_), sorted_ (true),
    04069             index1_data_ (capacity_), index2_data_ (capacity_), value_data_ (capacity_) {
    04070             storage_invariants ();
    04071             matrix_assign<scalar_assign> (*this, ae);
    04072         }
    04073
    04074         // Accessors
    04075         BOOST_UBLAS_INLINE
    04076         size_type size1 () const {
    04077             return size1_;
    04078         }
    04079         BOOST_UBLAS_INLINE
    04080         size_type size2 () const {
    04081             return size2_;
    04082         }
    04083         BOOST_UBLAS_INLINE
    04084         size_type nnz_capacity () const {
    04085             return capacity_;
    04086         }
    04087         BOOST_UBLAS_INLINE
    04088         size_type nnz () const {
    04089             return filled_;
    04090         }
    04091
    04092         // Storage accessors
    04093         BOOST_UBLAS_INLINE
    04094         static size_type index_base () {
    04095             return IB;
    04096         }
    04097         BOOST_UBLAS_INLINE
    04098         array_size_type filled () const {
    04099             return filled_;
    04100         }
    04101         BOOST_UBLAS_INLINE
    04102         const index_array_type &index1_data () const {
    04103             return index1_data_;
    04104         }
    04105         BOOST_UBLAS_INLINE
    04106         const index_array_type &index2_data () const {
    04107             return index2_data_;
    04108         }
    04109         BOOST_UBLAS_INLINE
    04110         const value_array_type &value_data () const {
    04111             return value_data_;
    04112         }
    04113         BOOST_UBLAS_INLINE
    04114         void set_filled (const array_size_type &filled) {
    04115             // Make sure that storage_invariants() succeeds
    04116             if (sorted_ && filled < filled_)
    04117                 sorted_filled_ = filled;
    04118             else
    04119                 sorted_ = (sorted_filled_ == filled);
    04120             filled_ = filled;
    04121             storage_invariants ();
    04122         }
    04123         BOOST_UBLAS_INLINE
    04124         index_array_type &index1_data () {
    04125             return index1_data_;
    04126         }
    04127         BOOST_UBLAS_INLINE
    04128         index_array_type &index2_data () {
    04129             return index2_data_;
    04130         }
    04131         BOOST_UBLAS_INLINE
    04132         value_array_type &value_data () {
    04133             return value_data_;
    04134         }
    04135
    04136         // Resizing
    04137     private:
    04138         BOOST_UBLAS_INLINE
    04139         array_size_type restrict_capacity (array_size_type non_zeros) const {
    04140             // minimum non_zeros
    04141             non_zeros = (std::max) (non_zeros, array_size_type((std::min) (size1_, size2_)));
    04142             // ISSUE no maximum as coordinate may contain inserted duplicates
    04143             return non_zeros;
    04144         }
    04145     public:
    04146         BOOST_UBLAS_INLINE
    04147         void resize (size_type size1, size_type size2, bool preserve = true) {
    04148             // FIXME preserve unimplemented
    04149             BOOST_UBLAS_CHECK (!preserve, internal_logic ());
    04150             size1_ = size1;
    04151             size2_ = size2;
    04152             capacity_ = restrict_capacity (capacity_);
    04153             index1_data_.resize (capacity_);
    04154             index2_data_.resize (capacity_);
    04155             value_data_.resize (capacity_);
    04156             filled_ = 0;
    04157             sorted_filled_ = filled_;
    04158             sorted_ = true;
    04159             storage_invariants ();
    04160         }
    04161
    04162         // Reserving
    04163         BOOST_UBLAS_INLINE
    04164         void reserve (array_size_type non_zeros, bool preserve = true) {
    04165             sort ();    // remove duplicate elements
    04166             capacity_ = restrict_capacity (non_zeros);
    04167             if (preserve) {
    04168                 index1_data_.resize (capacity_, size_type ());
    04169                 index2_data_.resize (capacity_, size_type ());
    04170                 value_data_.resize (capacity_, value_type ());
    04171                 filled_ = (std::min) (capacity_, filled_);
    04172             }
    04173             else {
    04174                 index1_data_.resize (capacity_);
    04175                 index2_data_.resize (capacity_);
    04176                 value_data_.resize (capacity_);
    04177                 filled_ = 0;
    04178             }
    04179             sorted_filled_ = filled_;
    04180             storage_invariants ();
    04181         }
    04182
    04183         // Element support
    04184         BOOST_UBLAS_INLINE
    04185         pointer find_element (size_type i, size_type j) {
    04186             return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i, j));
    04187         }
    04188         BOOST_UBLAS_INLINE
    04189         const_pointer find_element (size_type i, size_type j) const {
    04190             sort ();
    04191             size_type element1 (layout_type::index_M (i, j));
    04192             size_type element2 (layout_type::index_m (i, j));
    04193             vector_const_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (element1), std::less<size_type> ()));
    04194             vector_const_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (element1), std::less<size_type> ()));
    04195             if (itv_begin == itv_end)
    04196                 return 0;
    04197             const_subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ()));
    04198             const_subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ()));
    04199             const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (element2), std::less<size_type> ()));
    04200             if (it == it_end || *it != k_based (element2))
    04201                 return 0;
    04202             return &value_data_ [it - index2_data_.begin ()];
    04203         }
    04204
    04205         // Element access
    04206         BOOST_UBLAS_INLINE
    04207         const_reference operator () (size_type i, size_type j) const {
    04208             const_pointer p = find_element (i, j);
    04209             if (p)
    04210                 return *p;
    04211             else
    04212                 return zero_;
    04213         }
    04214         BOOST_UBLAS_INLINE
    04215         reference operator () (size_type i, size_type j) {
    04216 #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
    04217             pointer p = find_element (i, j);
    04218             if (p)
    04219                 return *p;
    04220             else
    04221                 return insert_element (i, j, value_type/*zero*/());
    04222 #else
    04223             return reference (*this, i, j);
    04224 #endif
    04225         }
    04226
    04227         // Element assignment
    04228         BOOST_UBLAS_INLINE
    04229         void append_element (size_type i, size_type j, const_reference t) {
    04230             if (filled_ >= capacity_)
    04231                 reserve (2 * filled_, true);
    04232             BOOST_UBLAS_CHECK (filled_ < capacity_, internal_logic ());
    04233             size_type element1 = layout_type::index_M (i, j);
    04234             size_type element2 = layout_type::index_m (i, j);
    04235             index1_data_ [filled_] = k_based (element1);
    04236             index2_data_ [filled_] = k_based (element2);
    04237             value_data_ [filled_] = t;
    04238             ++ filled_;
    04239             sorted_ = false;
    04240             storage_invariants ();
    04241         }
    04242         BOOST_UBLAS_INLINE
    04243         true_reference insert_element (size_type i, size_type j, const_reference t) {
    04244             BOOST_UBLAS_CHECK (!find_element (i, j), bad_index ());        // duplicate element
    04245             append_element (i, j, t);
    04246             return value_data_ [filled_ - 1];
    04247         }
    04248         BOOST_UBLAS_INLINE
    04249         void erase_element (size_type i, size_type j) {
    04250             size_type element1 = layout_type::index_M (i, j);
    04251             size_type element2 = layout_type::index_m (i, j);
    04252             sort ();
    04253             vector_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (element1), std::less<size_type> ()));
    04254             vector_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (element1), std::less<size_type> ()));
    04255             subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ()));
    04256             subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ()));
    04257             subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (element2), std::less<size_type> ()));
    04258             if (it != it_end && *it == k_based (element2)) {
    04259                 typename std::iterator_traits<subiterator_type>::difference_type n = it - index2_data_.begin ();
    04260                 vector_subiterator_type itv (index1_data_.begin () + n);
    04261                 std::copy (itv + 1, index1_data_.begin () + filled_, itv);
    04262                 std::copy (it + 1, index2_data_.begin () + filled_, it);
    04263                 typename value_array_type::iterator itt (value_data_.begin () + n);
    04264                 std::copy (itt + 1, value_data_.begin () + filled_, itt);
    04265                 -- filled_;
    04266                 sorted_filled_ = filled_;
    04267             }
    04268             storage_invariants ();
    04269         }
    04270
    04271         // Zeroing
    04272         BOOST_UBLAS_INLINE
    04273         void clear () {
    04274             filled_ = 0;
    04275             sorted_filled_ = filled_;
    04276             sorted_ = true;
    04277             storage_invariants ();
    04278         }
    04279
    04280         // Assignment
    04281         BOOST_UBLAS_INLINE
    04282         coordinate_matrix &operator = (const coordinate_matrix &m) {
    04283             if (this != &m) {
    04284                 size1_ = m.size1_;
    04285                 size2_ = m.size2_;
    04286                 capacity_ = m.capacity_;
    04287                 filled_ = m.filled_;
    04288                 sorted_filled_ = m.sorted_filled_;
    04289                 sorted_ = m.sorted_;
    04290                 index1_data_ = m.index1_data_;
    04291                 index2_data_ = m.index2_data_;
    04292                 value_data_ = m.value_data_;
    04293                 BOOST_UBLAS_CHECK (capacity_ == index1_data_.size (), internal_logic ());
    04294                 BOOST_UBLAS_CHECK (capacity_ == index2_data_.size (), internal_logic ());
    04295                 BOOST_UBLAS_CHECK (capacity_ == value_data_.size (), internal_logic ());
    04296             }
    04297             storage_invariants ();
    04298             return *this;
    04299         }
    04300         template<class C>          // Container assignment without temporary
    04301         BOOST_UBLAS_INLINE
    04302         coordinate_matrix &operator = (const matrix_container<C> &m) {
    04303             resize (m ().size1 (), m ().size2 (), false);
    04304             assign (m);
    04305             return *this;
    04306         }
    04307         BOOST_UBLAS_INLINE
    04308         coordinate_matrix &assign_temporary (coordinate_matrix &m) {
    04309             swap (m);
    04310             return *this;
    04311         }
    04312         template<class AE>
    04313         BOOST_UBLAS_INLINE
    04314         coordinate_matrix &operator = (const matrix_expression<AE> &ae) {
    04315             self_type temporary (ae, capacity_);
    04316             return assign_temporary (temporary);
    04317         }
    04318         template<class AE>
    04319         BOOST_UBLAS_INLINE
    04320         coordinate_matrix &assign (const matrix_expression<AE> &ae) {
    04321             matrix_assign<scalar_assign> (*this, ae);
    04322             return *this;
    04323         }
    04324         template<class AE>
    04325         BOOST_UBLAS_INLINE
    04326         coordinate_matrix& operator += (const matrix_expression<AE> &ae) {
    04327             self_type temporary (*this + ae, capacity_);
    04328             return assign_temporary (temporary);
    04329         }
    04330         template<class C>          // Container assignment without temporary
    04331         BOOST_UBLAS_INLINE
    04332         coordinate_matrix &operator += (const matrix_container<C> &m) {
    04333             plus_assign (m);
    04334             return *this;
    04335         }
    04336         template<class AE>
    04337         BOOST_UBLAS_INLINE
    04338         coordinate_matrix &plus_assign (const matrix_expression<AE> &ae) {
    04339             matrix_assign<scalar_plus_assign> (*this, ae);
    04340             return *this;
    04341         }
    04342         template<class AE>
    04343         BOOST_UBLAS_INLINE
    04344         coordinate_matrix& operator -= (const matrix_expression<AE> &ae) {
    04345             self_type temporary (*this - ae, capacity_);
    04346             return assign_temporary (temporary);
    04347         }
    04348         template<class C>          // Container assignment without temporary
    04349         BOOST_UBLAS_INLINE
    04350         coordinate_matrix &operator -= (const matrix_container<C> &m) {
    04351             minus_assign (m);
    04352             return *this;
    04353         }
    04354         template<class AE>
    04355         BOOST_UBLAS_INLINE
    04356         coordinate_matrix &minus_assign (const matrix_expression<AE> &ae) {
    04357             matrix_assign<scalar_minus_assign> (*this, ae);
    04358             return *this;
    04359         }
    04360         template<class AT>
    04361         BOOST_UBLAS_INLINE
    04362         coordinate_matrix& operator *= (const AT &at) {
    04363             matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
    04364             return *this;
    04365         }
    04366         template<class AT>
    04367         BOOST_UBLAS_INLINE
    04368         coordinate_matrix& operator /= (const AT &at) {
    04369             matrix_assign_scalar<scalar_divides_assign> (*this, at);
    04370             return *this;
    04371         }
    04372
    04373         // Swapping
    04374         BOOST_UBLAS_INLINE
    04375         void swap (coordinate_matrix &m) {
    04376             if (this != &m) {
    04377                 std::swap (size1_, m.size1_);
    04378                 std::swap (size2_, m.size2_);
    04379                 std::swap (capacity_, m.capacity_);
    04380                 std::swap (filled_, m.filled_);
    04381                 std::swap (sorted_filled_, m.sorted_filled_);
    04382                 std::swap (sorted_, m.sorted_);
    04383                 index1_data_.swap (m.index1_data_);
    04384                 index2_data_.swap (m.index2_data_);
    04385                 value_data_.swap (m.value_data_);
    04386             }
    04387             storage_invariants ();
    04388         }
    04389         BOOST_UBLAS_INLINE
    04390         friend void swap (coordinate_matrix &m1, coordinate_matrix &m2) {
    04391             m1.swap (m2);
    04392         }
    04393
    04394         // Sorting and summation of duplicates
    04395         BOOST_UBLAS_INLINE
    04396         void sort () const {
    04397             if (! sorted_ && filled_ > 0) {
    04398                 typedef index_triple_array<index_array_type, index_array_type, value_array_type> array_triple;
    04399                 array_triple ita (filled_, index1_data_, index2_data_, value_data_);
    04400                 const typename array_triple::iterator iunsorted = ita.begin () + sorted_filled_;
    04401                 // sort new elements and merge
    04402                 std::sort (iunsorted, ita.end ());
    04403                 std::inplace_merge (ita.begin (), iunsorted, ita.end ());
    04404
    04405                 // sum duplicates with += and remove
    04406                 array_size_type filled = 0;
    04407                 for (array_size_type i = 1; i < filled_; ++ i) {
    04408                     if (index1_data_ [filled] != index1_data_ [i] ||
    04409                         index2_data_ [filled] != index2_data_ [i]) {
    04410                         ++ filled;
    04411                         if (filled != i) {
    04412                             index1_data_ [filled] = index1_data_ [i];
    04413                             index2_data_ [filled] = index2_data_ [i];
    04414                             value_data_ [filled] = value_data_ [i];
    04415                         }
    04416                     } else {
    04417                         value_data_ [filled] += value_data_ [i];
    04418                     }
    04419                 }
    04420                 filled_ = filled + 1;
    04421                 sorted_filled_ = filled_;
    04422                 sorted_ = true;
    04423                 storage_invariants ();
    04424             }
    04425         }
    04426
    04427         // Back element insertion and erasure
    04428         BOOST_UBLAS_INLINE
    04429         void push_back (size_type i, size_type j, const_reference t) {
    04430             size_type element1 = layout_type::index_M (i, j);
    04431             size_type element2 = layout_type::index_m (i, j);
    04432             // must maintain sort order
    04433             BOOST_UBLAS_CHECK (sorted_ &&
    04434                     (filled_ == 0 ||
    04435                     index1_data_ [filled_ - 1] < k_based (element1) ||
    04436                     (index1_data_ [filled_ - 1] == k_based (element1) && index2_data_ [filled_ - 1] < k_based (element2)))
    04437                     , external_logic ());
    04438             if (filled_ >= capacity_)
    04439                 reserve (2 * filled_, true);
    04440             BOOST_UBLAS_CHECK (filled_ < capacity_, internal_logic ());
    04441             index1_data_ [filled_] = k_based (element1);
    04442             index2_data_ [filled_] = k_based (element2);
    04443             value_data_ [filled_] = t;
    04444             ++ filled_;
    04445             sorted_filled_ = filled_;
    04446             storage_invariants ();
    04447         }
    04448         BOOST_UBLAS_INLINE
    04449         void pop_back () {
    04450             // ISSUE invariants could be simpilfied if sorted required as precondition
    04451             BOOST_UBLAS_CHECK (filled_ > 0, external_logic ());
    04452             -- filled_;
    04453             sorted_filled_ = (std::min) (sorted_filled_, filled_);
    04454             sorted_ = sorted_filled_ = filled_;
    04455             storage_invariants ();
    04456         }
    04457
    04458         // Iterator types
    04459     private:
    04460         // Use index array iterator
    04461         typedef typename IA::const_iterator vector_const_subiterator_type;
    04462         typedef typename IA::iterator vector_subiterator_type;
    04463         typedef typename IA::const_iterator const_subiterator_type;
    04464         typedef typename IA::iterator subiterator_type;
    04465
    04466         BOOST_UBLAS_INLINE
    04467         true_reference at_element (size_type i, size_type j) {
    04468             pointer p = find_element (i, j);
    04469             BOOST_UBLAS_CHECK (p, bad_index ());
    04470             return *p;
    04471         }
    04472
    04473     public:
    04474         class const_iterator1;
    04475         class iterator1;
    04476         class const_iterator2;
    04477         class iterator2;
    04478         typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
    04479         typedef reverse_iterator_base1<iterator1> reverse_iterator1;
    04480         typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
    04481         typedef reverse_iterator_base2<iterator2> reverse_iterator2;
    04482
    04483         // Element lookup
    04484         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    04485         const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const {
    04486             sort ();
    04487             for (;;) {
    04488                 size_type address1 (layout_type::index_M (i, j));
    04489                 size_type address2 (layout_type::index_m (i, j));
    04490                 vector_const_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ()));
    04491                 vector_const_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ()));
    04492
    04493                 const_subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ()));
    04494                 const_subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ()));
    04495
    04496                 const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ()));
    04497                 vector_const_subiterator_type itv (index1_data_.begin () + (it - index2_data_.begin ()));
    04498                 if (rank == 0)
    04499                     return const_iterator1 (*this, rank, i, j, itv, it);
    04500                 if (it != it_end && zero_based (*it) == address2)
    04501                     return const_iterator1 (*this, rank, i, j, itv, it);
    04502                 if (direction > 0) {
    04503                     if (layout_type::fast_i ()) {
    04504                         if (it == it_end)
    04505                             return const_iterator1 (*this, rank, i, j, itv, it);
    04506                         i = zero_based (*it);
    04507                     } else {
    04508                         if (i >= size1_)
    04509                             return const_iterator1 (*this, rank, i, j, itv, it);
    04510                         ++ i;
    04511                     }
    04512                 } else /* if (direction < 0)  */ {
    04513                     if (layout_type::fast_i ()) {
    04514                         if (it == index2_data_.begin () + array_size_type (zero_based (*itv)))
    04515                             return const_iterator1 (*this, rank, i, j, itv, it);
    04516                         i = zero_based (*(it - 1));
    04517                     } else {
    04518                         if (i == 0)
    04519                             return const_iterator1 (*this, rank, i, j, itv, it);
    04520                         -- i;
    04521                     }
    04522                 }
    04523             }
    04524         }
    04525         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    04526         iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) {
    04527             sort ();
    04528             for (;;) {
    04529                 size_type address1 (layout_type::index_M (i, j));
    04530                 size_type address2 (layout_type::index_m (i, j));
    04531                 vector_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ()));
    04532                 vector_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ()));
    04533
    04534                 subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ()));
    04535                 subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ()));
    04536
    04537                 subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ()));
    04538                 vector_subiterator_type itv (index1_data_.begin () + (it - index2_data_.begin ()));
    04539                 if (rank == 0)
    04540                     return iterator1 (*this, rank, i, j, itv, it);
    04541                 if (it != it_end && zero_based (*it) == address2)
    04542                     return iterator1 (*this, rank, i, j, itv, it);
    04543                 if (direction > 0) {
    04544                     if (layout_type::fast_i ()) {
    04545                         if (it == it_end)
    04546                             return iterator1 (*this, rank, i, j, itv, it);
    04547                         i = zero_based (*it);
    04548                     } else {
    04549                         if (i >= size1_)
    04550                             return iterator1 (*this, rank, i, j, itv, it);
    04551                         ++ i;
    04552                     }
    04553                 } else /* if (direction < 0)  */ {
    04554                     if (layout_type::fast_i ()) {
    04555                         if (it == index2_data_.begin () + array_size_type (zero_based (*itv)))
    04556                             return iterator1 (*this, rank, i, j, itv, it);
    04557                         i = zero_based (*(it - 1));
    04558                     } else {
    04559                         if (i == 0)
    04560                             return iterator1 (*this, rank, i, j, itv, it);
    04561                         -- i;
    04562                     }
    04563                 }
    04564             }
    04565         }
    04566         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    04567         const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const {
    04568             sort ();
    04569             for (;;) {
    04570                 size_type address1 (layout_type::index_M (i, j));
    04571                 size_type address2 (layout_type::index_m (i, j));
    04572                 vector_const_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ()));
    04573                 vector_const_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ()));
    04574
    04575                 const_subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ()));
    04576                 const_subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ()));
    04577
    04578                 const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ()));
    04579                 vector_const_subiterator_type itv (index1_data_.begin () + (it - index2_data_.begin ()));
    04580                 if (rank == 0)
    04581                     return const_iterator2 (*this, rank, i, j, itv, it);
    04582                 if (it != it_end && zero_based (*it) == address2)
    04583                     return const_iterator2 (*this, rank, i, j, itv, it);
    04584                 if (direction > 0) {
    04585                     if (layout_type::fast_j ()) {
    04586                         if (it == it_end)
    04587                             return const_iterator2 (*this, rank, i, j, itv, it);
    04588                         j = zero_based (*it);
    04589                     } else {
    04590                         if (j >= size2_)
    04591                             return const_iterator2 (*this, rank, i, j, itv, it);
    04592                         ++ j;
    04593                     }
    04594                 } else /* if (direction < 0)  */ {
    04595                     if (layout_type::fast_j ()) {
    04596                         if (it == index2_data_.begin () + array_size_type (zero_based (*itv)))
    04597                             return const_iterator2 (*this, rank, i, j, itv, it);
    04598                         j = zero_based (*(it - 1));
    04599                     } else {
    04600                         if (j == 0)
    04601                             return const_iterator2 (*this, rank, i, j, itv, it);
    04602                         -- j;
    04603                     }
    04604                 }
    04605             }
    04606         }
    04607         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    04608         iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) {
    04609             sort ();
    04610             for (;;) {
    04611                 size_type address1 (layout_type::index_M (i, j));
    04612                 size_type address2 (layout_type::index_m (i, j));
    04613                 vector_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ()));
    04614                 vector_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ()));
    04615
    04616                 subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ()));
    04617                 subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ()));
    04618
    04619                 subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ()));
    04620                 vector_subiterator_type itv (index1_data_.begin () + (it - index2_data_.begin ()));
    04621                 if (rank == 0)
    04622                     return iterator2 (*this, rank, i, j, itv, it);
    04623                 if (it != it_end && zero_based (*it) == address2)
    04624                     return iterator2 (*this, rank, i, j, itv, it);
    04625                 if (direction > 0) {
    04626                     if (layout_type::fast_j ()) {
    04627                         if (it == it_end)
    04628                             return iterator2 (*this, rank, i, j, itv, it);
    04629                         j = zero_based (*it);
    04630                     } else {
    04631                         if (j >= size2_)
    04632                             return iterator2 (*this, rank, i, j, itv, it);
    04633                         ++ j;
    04634                     }
    04635                 } else /* if (direction < 0)  */ {
    04636                     if (layout_type::fast_j ()) {
    04637                         if (it == index2_data_.begin () + array_size_type (zero_based (*itv)))
    04638                             return iterator2 (*this, rank, i, j, itv, it);
    04639                         j = zero_based (*(it - 1));
    04640                     } else {
    04641                         if (j == 0)
    04642                             return iterator2 (*this, rank, i, j, itv, it);
    04643                         -- j;
    04644                     }
    04645                 }
    04646             }
    04647         }
    04648
    04649
    04650         class const_iterator1:
    04651             public container_const_reference<coordinate_matrix>,
    04652             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
    04653                                                const_iterator1, value_type> {
    04654         public:
    04655             typedef typename coordinate_matrix::value_type value_type;
    04656             typedef typename coordinate_matrix::difference_type difference_type;
    04657             typedef typename coordinate_matrix::const_reference reference;
    04658             typedef const typename coordinate_matrix::pointer pointer;
    04659
    04660             typedef const_iterator2 dual_iterator_type;
    04661             typedef const_reverse_iterator2 dual_reverse_iterator_type;
    04662
    04663             // Construction and destruction
    04664             BOOST_UBLAS_INLINE
    04665             const_iterator1 ():
    04666                 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
    04667             BOOST_UBLAS_INLINE
    04668             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):
    04669                 container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
    04670             BOOST_UBLAS_INLINE
    04671             const_iterator1 (const iterator1 &it):
    04672                 container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {}
    04673
    04674             // Arithmetic
    04675             BOOST_UBLAS_INLINE
    04676             const_iterator1 &operator ++ () {
    04677                 if (rank_ == 1 && layout_type::fast_i ())
    04678                     ++ it_;
    04679                 else {
    04680                     i_ = index1 () + 1;
    04681                     if (rank_ == 1)
    04682                         *this = (*this) ().find1 (rank_, i_, j_, 1);
    04683                 }
    04684                 return *this;
    04685             }
    04686             BOOST_UBLAS_INLINE
    04687             const_iterator1 &operator -- () {
    04688                 if (rank_ == 1 && layout_type::fast_i ())
    04689                     -- it_;
    04690                 else {
    04691                     i_ = index1 () - 1;
    04692                     if (rank_ == 1)
    04693                         *this = (*this) ().find1 (rank_, i_, j_, -1);
    04694                 }
    04695                 return *this;
    04696             }
    04697
    04698             // Dereference
    04699             BOOST_UBLAS_INLINE
    04700             const_reference operator * () const {
    04701                 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
    04702                 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
    04703                 if (rank_ == 1) {
    04704                     return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
    04705                 } else {
    04706                     return (*this) () (i_, j_);
    04707                 }
    04708             }
    04709
    04710 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    04711             BOOST_UBLAS_INLINE
    04712 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    04713             typename self_type::
    04714 #endif
    04715             const_iterator2 begin () const {
    04716                 const self_type &m = (*this) ();
    04717                 return m.find2 (1, index1 (), 0);
    04718             }
    04719             BOOST_UBLAS_INLINE
    04720 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    04721             typename self_type::
    04722 #endif
    04723             const_iterator2 end () const {
    04724                 const self_type &m = (*this) ();
    04725                 return m.find2 (1, index1 (), m.size2 ());
    04726             }
    04727             BOOST_UBLAS_INLINE
    04728 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    04729             typename self_type::
    04730 #endif
    04731             const_reverse_iterator2 rbegin () const {
    04732                 return const_reverse_iterator2 (end ());
    04733             }
    04734             BOOST_UBLAS_INLINE
    04735 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    04736             typename self_type::
    04737 #endif
    04738             const_reverse_iterator2 rend () const {
    04739                 return const_reverse_iterator2 (begin ());
    04740             }
    04741 #endif
    04742
    04743             // Indices
    04744             BOOST_UBLAS_INLINE
    04745             size_type index1 () const {
    04746                 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
    04747                 if (rank_ == 1) {
    04748                     BOOST_UBLAS_CHECK (layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
    04749                     return layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
    04750                 } else {
    04751                     return i_;
    04752                 }
    04753             }
    04754             BOOST_UBLAS_INLINE
    04755             size_type index2 () const {
    04756                 if (rank_ == 1) {
    04757                     BOOST_UBLAS_CHECK (layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
    04758                     return layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
    04759                 } else {
    04760                     return j_;
    04761                 }
    04762             }
    04763
    04764             // Assignment
    04765             BOOST_UBLAS_INLINE
    04766             const_iterator1 &operator = (const const_iterator1 &it) {
    04767                 container_const_reference<self_type>::assign (&it ());
    04768                 rank_ = it.rank_;
    04769                 i_ = it.i_;
    04770                 j_ = it.j_;
    04771                 itv_ = it.itv_;
    04772                 it_ = it.it_;
    04773                 return *this;
    04774             }
    04775
    04776             // Comparison
    04777             BOOST_UBLAS_INLINE
    04778             bool operator == (const const_iterator1 &it) const {
    04779                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
    04780                 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
    04781                 if (rank_ == 1 || it.rank_ == 1) {
    04782                     return it_ == it.it_;
    04783                 } else {
    04784                     return i_ == it.i_ && j_ == it.j_;
    04785                 }
    04786             }
    04787
    04788         private:
    04789             int rank_;
    04790             size_type i_;
    04791             size_type j_;
    04792             vector_const_subiterator_type itv_;
    04793             const_subiterator_type it_;
    04794         };
    04795
    04796         BOOST_UBLAS_INLINE
    04797         const_iterator1 begin1 () const {
    04798             return find1 (0, 0, 0);
    04799         }
    04800         BOOST_UBLAS_INLINE
    04801         const_iterator1 end1 () const {
    04802             return find1 (0, size1_, 0);
    04803         }
    04804
    04805         class iterator1:
    04806             public container_reference<coordinate_matrix>,
    04807             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
    04808                                                iterator1, value_type> {
    04809         public:
    04810             typedef typename coordinate_matrix::value_type value_type;
    04811             typedef typename coordinate_matrix::difference_type difference_type;
    04812             typedef typename coordinate_matrix::true_reference reference;
    04813             typedef typename coordinate_matrix::pointer pointer;
    04814
    04815             typedef iterator2 dual_iterator_type;
    04816             typedef reverse_iterator2 dual_reverse_iterator_type;
    04817
    04818             // Construction and destruction
    04819             BOOST_UBLAS_INLINE
    04820             iterator1 ():
    04821                 container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
    04822             BOOST_UBLAS_INLINE
    04823             iterator1 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it):
    04824                 container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
    04825
    04826             // Arithmetic
    04827             BOOST_UBLAS_INLINE
    04828             iterator1 &operator ++ () {
    04829                 if (rank_ == 1 && layout_type::fast_i ())
    04830                     ++ it_;
    04831                 else {
    04832                     i_ = index1 () + 1;
    04833                     if (rank_ == 1)
    04834                         *this = (*this) ().find1 (rank_, i_, j_, 1);
    04835                 }
    04836                 return *this;
    04837             }
    04838             BOOST_UBLAS_INLINE
    04839             iterator1 &operator -- () {
    04840                 if (rank_ == 1 && layout_type::fast_i ())
    04841                     -- it_;
    04842                 else {
    04843                     i_ = index1 () - 1;
    04844                     if (rank_ == 1)
    04845                         *this = (*this) ().find1 (rank_, i_, j_, -1);
    04846                 }
    04847                 return *this;
    04848             }
    04849
    04850             // Dereference
    04851             BOOST_UBLAS_INLINE
    04852             reference operator * () const {
    04853                 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
    04854                 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
    04855                 if (rank_ == 1) {
    04856                     return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
    04857                 } else {
    04858                     return (*this) ().at_element (i_, j_);
    04859                 }
    04860             }
    04861
    04862 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    04863             BOOST_UBLAS_INLINE
    04864 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    04865             typename self_type::
    04866 #endif
    04867             iterator2 begin () const {
    04868                 self_type &m = (*this) ();
    04869                 return m.find2 (1, index1 (), 0);
    04870             }
    04871             BOOST_UBLAS_INLINE
    04872 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    04873             typename self_type::
    04874 #endif
    04875             iterator2 end () const {
    04876                 self_type &m = (*this) ();
    04877                 return m.find2 (1, index1 (), m.size2 ());
    04878             }
    04879             BOOST_UBLAS_INLINE
    04880 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    04881             typename self_type::
    04882 #endif
    04883             reverse_iterator2 rbegin () const {
    04884                 return reverse_iterator2 (end ());
    04885             }
    04886             BOOST_UBLAS_INLINE
    04887 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    04888             typename self_type::
    04889 #endif
    04890             reverse_iterator2 rend () const {
    04891                 return reverse_iterator2 (begin ());
    04892             }
    04893 #endif
    04894
    04895             // Indices
    04896             BOOST_UBLAS_INLINE
    04897             size_type index1 () const {
    04898                 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
    04899                 if (rank_ == 1) {
    04900                     BOOST_UBLAS_CHECK (layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
    04901                     return layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
    04902                 } else {
    04903                     return i_;
    04904                 }
    04905             }
    04906             BOOST_UBLAS_INLINE
    04907             size_type index2 () const {
    04908                 if (rank_ == 1) {
    04909                     BOOST_UBLAS_CHECK (layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
    04910                     return layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
    04911                 } else {
    04912                     return j_;
    04913                 }
    04914             }
    04915
    04916             // Assignment
    04917             BOOST_UBLAS_INLINE
    04918             iterator1 &operator = (const iterator1 &it) {
    04919                 container_reference<self_type>::assign (&it ());
    04920                 rank_ = it.rank_;
    04921                 i_ = it.i_;
    04922                 j_ = it.j_;
    04923                 itv_ = it.itv_;
    04924                 it_ = it.it_;
    04925                 return *this;
    04926             }
    04927
    04928             // Comparison
    04929             BOOST_UBLAS_INLINE
    04930             bool operator == (const iterator1 &it) const {
    04931                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
    04932                 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
    04933                 if (rank_ == 1 || it.rank_ == 1) {
    04934                     return it_ == it.it_;
    04935                 } else {
    04936                     return i_ == it.i_ && j_ == it.j_;
    04937                 }
    04938             }
    04939
    04940         private:
    04941             int rank_;
    04942             size_type i_;
    04943             size_type j_;
    04944             vector_subiterator_type itv_;
    04945             subiterator_type it_;
    04946
    04947             friend class const_iterator1;
    04948         };
    04949
    04950         BOOST_UBLAS_INLINE
    04951         iterator1 begin1 () {
    04952             return find1 (0, 0, 0);
    04953         }
    04954         BOOST_UBLAS_INLINE
    04955         iterator1 end1 () {
    04956             return find1 (0, size1_, 0);
    04957         }
    04958
    04959         class const_iterator2:
    04960             public container_const_reference<coordinate_matrix>,
    04961             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
    04962                                                const_iterator2, value_type> {
    04963         public:
    04964             typedef typename coordinate_matrix::value_type value_type;
    04965             typedef typename coordinate_matrix::difference_type difference_type;
    04966             typedef typename coordinate_matrix::const_reference reference;
    04967             typedef const typename coordinate_matrix::pointer pointer;
    04968
    04969             typedef const_iterator1 dual_iterator_type;
    04970             typedef const_reverse_iterator1 dual_reverse_iterator_type;
    04971
    04972             // Construction and destruction
    04973             BOOST_UBLAS_INLINE
    04974             const_iterator2 ():
    04975                 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
    04976             BOOST_UBLAS_INLINE
    04977             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):
    04978                 container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
    04979             BOOST_UBLAS_INLINE
    04980             const_iterator2 (const iterator2 &it):
    04981                 container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {}
    04982
    04983             // Arithmetic
    04984             BOOST_UBLAS_INLINE
    04985             const_iterator2 &operator ++ () {
    04986                 if (rank_ == 1 && layout_type::fast_j ())
    04987                     ++ it_;
    04988                 else {
    04989                     j_ = index2 () + 1;
    04990                     if (rank_ == 1)
    04991                         *this = (*this) ().find2 (rank_, i_, j_, 1);
    04992                 }
    04993                 return *this;
    04994             }
    04995             BOOST_UBLAS_INLINE
    04996             const_iterator2 &operator -- () {
    04997                 if (rank_ == 1 && layout_type::fast_j ())
    04998                     -- it_;
    04999                 else {
    05000                     j_ = index2 () - 1;
    05001                     if (rank_ == 1)
    05002                         *this = (*this) ().find2 (rank_, i_, j_, -1);
    05003                 }
    05004                 return *this;
    05005             }
    05006
    05007             // Dereference
    05008             BOOST_UBLAS_INLINE
    05009             const_reference operator * () const {
    05010                 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
    05011                 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
    05012                 if (rank_ == 1) {
    05013                     return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
    05014                 } else {
    05015                     return (*this) () (i_, j_);
    05016                 }
    05017             }
    05018
    05019 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    05020             BOOST_UBLAS_INLINE
    05021 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    05022             typename self_type::
    05023 #endif
    05024             const_iterator1 begin () const {
    05025                 const self_type &m = (*this) ();
    05026                 return m.find1 (1, 0, index2 ());
    05027             }
    05028             BOOST_UBLAS_INLINE
    05029 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    05030             typename self_type::
    05031 #endif
    05032             const_iterator1 end () const {
    05033                 const self_type &m = (*this) ();
    05034                 return m.find1 (1, m.size1 (), index2 ());
    05035             }
    05036             BOOST_UBLAS_INLINE
    05037 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    05038             typename self_type::
    05039 #endif
    05040             const_reverse_iterator1 rbegin () const {
    05041                 return const_reverse_iterator1 (end ());
    05042             }
    05043             BOOST_UBLAS_INLINE
    05044 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    05045             typename self_type::
    05046 #endif
    05047             const_reverse_iterator1 rend () const {
    05048                 return const_reverse_iterator1 (begin ());
    05049             }
    05050 #endif
    05051
    05052             // Indices
    05053             BOOST_UBLAS_INLINE
    05054             size_type index1 () const {
    05055                 if (rank_ == 1) {
    05056                     BOOST_UBLAS_CHECK (layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
    05057                     return layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
    05058                 } else {
    05059                     return i_;
    05060                 }
    05061             }
    05062             BOOST_UBLAS_INLINE
    05063             size_type index2 () const {
    05064                 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
    05065                 if (rank_ == 1) {
    05066                     BOOST_UBLAS_CHECK (layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
    05067                     return layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
    05068                 } else {
    05069                     return j_;
    05070                 }
    05071             }
    05072
    05073             // Assignment
    05074             BOOST_UBLAS_INLINE
    05075             const_iterator2 &operator = (const const_iterator2 &it) {
    05076                 container_const_reference<self_type>::assign (&it ());
    05077                 rank_ = it.rank_;
    05078                 i_ = it.i_;
    05079                 j_ = it.j_;
    05080                 itv_ = it.itv_;
    05081                 it_ = it.it_;
    05082                 return *this;
    05083             }
    05084
    05085             // Comparison
    05086             BOOST_UBLAS_INLINE
    05087             bool operator == (const const_iterator2 &it) const {
    05088                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
    05089                 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
    05090                 if (rank_ == 1 || it.rank_ == 1) {
    05091                     return it_ == it.it_;
    05092                 } else {
    05093                     return i_ == it.i_ && j_ == it.j_;
    05094                 }
    05095             }
    05096
    05097         private:
    05098             int rank_;
    05099             size_type i_;
    05100             size_type j_;
    05101             vector_const_subiterator_type itv_;
    05102             const_subiterator_type it_;
    05103         };
    05104
    05105         BOOST_UBLAS_INLINE
    05106         const_iterator2 begin2 () const {
    05107             return find2 (0, 0, 0);
    05108         }
    05109         BOOST_UBLAS_INLINE
    05110         const_iterator2 end2 () const {
    05111             return find2 (0, 0, size2_);
    05112         }
    05113
    05114         class iterator2:
    05115             public container_reference<coordinate_matrix>,
    05116             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
    05117                                                iterator2, value_type> {
    05118         public:
    05119             typedef typename coordinate_matrix::value_type value_type;
    05120             typedef typename coordinate_matrix::difference_type difference_type;
    05121             typedef typename coordinate_matrix::true_reference reference;
    05122             typedef typename coordinate_matrix::pointer pointer;
    05123
    05124             typedef iterator1 dual_iterator_type;
    05125             typedef reverse_iterator1 dual_reverse_iterator_type;
    05126
    05127             // Construction and destruction
    05128             BOOST_UBLAS_INLINE
    05129             iterator2 ():
    05130                 container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
    05131             BOOST_UBLAS_INLINE
    05132             iterator2 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it):
    05133                 container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
    05134
    05135             // Arithmetic
    05136             BOOST_UBLAS_INLINE
    05137             iterator2 &operator ++ () {
    05138                 if (rank_ == 1 && layout_type::fast_j ())
    05139                     ++ it_;
    05140                 else {
    05141                     j_ = index2 () + 1;
    05142                     if (rank_ == 1)
    05143                         *this = (*this) ().find2 (rank_, i_, j_, 1);
    05144                 }
    05145                 return *this;
    05146             }
    05147             BOOST_UBLAS_INLINE
    05148             iterator2 &operator -- () {
    05149                 if (rank_ == 1 && layout_type::fast_j ())
    05150                     -- it_;
    05151                 else {
    05152                     j_ = index2 ();
    05153                     if (rank_ == 1)
    05154                         *this = (*this) ().find2 (rank_, i_, j_, -1);
    05155                 }
    05156                 return *this;
    05157             }
    05158
    05159             // Dereference
    05160             BOOST_UBLAS_INLINE
    05161             reference operator * () const {
    05162                 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
    05163                 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
    05164                 if (rank_ == 1) {
    05165                     return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
    05166                 } else {
    05167                     return (*this) ().at_element (i_, j_);
    05168                 }
    05169             }
    05170
    05171 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
    05172             BOOST_UBLAS_INLINE
    05173 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    05174             typename self_type::
    05175 #endif
    05176             iterator1 begin () const {
    05177                 self_type &m = (*this) ();
    05178                 return m.find1 (1, 0, index2 ());
    05179             }
    05180             BOOST_UBLAS_INLINE
    05181 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    05182             typename self_type::
    05183 #endif
    05184             iterator1 end () const {
    05185                 self_type &m = (*this) ();
    05186                 return m.find1 (1, m.size1 (), index2 ());
    05187             }
    05188             BOOST_UBLAS_INLINE
    05189 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    05190             typename self_type::
    05191 #endif
    05192             reverse_iterator1 rbegin () const {
    05193                 return reverse_iterator1 (end ());
    05194             }
    05195             BOOST_UBLAS_INLINE
    05196 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
    05197             typename self_type::
    05198 #endif
    05199             reverse_iterator1 rend () const {
    05200                 return reverse_iterator1 (begin ());
    05201             }
    05202 #endif
    05203
    05204             // Indices
    05205             BOOST_UBLAS_INLINE
    05206             size_type index1 () const {
    05207                 if (rank_ == 1) {
    05208                     BOOST_UBLAS_CHECK (layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
    05209                     return layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
    05210                 } else {
    05211                     return i_;
    05212                 }
    05213             }
    05214             BOOST_UBLAS_INLINE
    05215             size_type index2 () const {
    05216                 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
    05217                 if (rank_ == 1) {
    05218                     BOOST_UBLAS_CHECK (layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
    05219                     return layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
    05220                 } else {
    05221                     return j_;
    05222                 }
    05223             }
    05224
    05225             // Assignment
    05226             BOOST_UBLAS_INLINE
    05227             iterator2 &operator = (const iterator2 &it) {
    05228                 container_reference<self_type>::assign (&it ());
    05229                 rank_ = it.rank_;
    05230                 i_ = it.i_;
    05231                 j_ = it.j_;
    05232                 itv_ = it.itv_;
    05233                 it_ = it.it_;
    05234                 return *this;
    05235             }
    05236
    05237             // Comparison
    05238             BOOST_UBLAS_INLINE
    05239             bool operator == (const iterator2 &it) const {
    05240                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
    05241                 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
    05242                 if (rank_ == 1 || it.rank_ == 1) {
    05243                     return it_ == it.it_;
    05244                 } else {
    05245                     return i_ == it.i_ && j_ == it.j_;
    05246                 }
    05247             }
    05248
    05249         private:
    05250             int rank_;
    05251             size_type i_;
    05252             size_type j_;
    05253             vector_subiterator_type itv_;
    05254             subiterator_type it_;
    05255
    05256             friend class const_iterator2;
    05257         };
    05258
    05259         BOOST_UBLAS_INLINE
    05260         iterator2 begin2 () {
    05261             return find2 (0, 0, 0);
    05262         }
    05263         BOOST_UBLAS_INLINE
    05264         iterator2 end2 () {
    05265             return find2 (0, 0, size2_);
    05266         }
    05267
    05268         // Reverse iterators
    05269
    05270         BOOST_UBLAS_INLINE
    05271         const_reverse_iterator1 rbegin1 () const {
    05272             return const_reverse_iterator1 (end1 ());
    05273         }
    05274         BOOST_UBLAS_INLINE
    05275         const_reverse_iterator1 rend1 () const {
    05276             return const_reverse_iterator1 (begin1 ());
    05277         }
    05278
    05279         BOOST_UBLAS_INLINE
    05280         reverse_iterator1 rbegin1 () {
    05281             return reverse_iterator1 (end1 ());
    05282         }
    05283         BOOST_UBLAS_INLINE
    05284         reverse_iterator1 rend1 () {
    05285             return reverse_iterator1 (begin1 ());
    05286         }
    05287
    05288         BOOST_UBLAS_INLINE
    05289         const_reverse_iterator2 rbegin2 () const {
    05290             return const_reverse_iterator2 (end2 ());
    05291         }
    05292         BOOST_UBLAS_INLINE
    05293         const_reverse_iterator2 rend2 () const {
    05294             return const_reverse_iterator2 (begin2 ());
    05295         }
    05296
    05297         BOOST_UBLAS_INLINE
    05298         reverse_iterator2 rbegin2 () {
    05299             return reverse_iterator2 (end2 ());
    05300         }
    05301         BOOST_UBLAS_INLINE
    05302         reverse_iterator2 rend2 () {
    05303             return reverse_iterator2 (begin2 ());
    05304         }
    05305
    05306          // Serialization
    05307         template<class Archive>
    05308         void serialize(Archive & ar, const unsigned int /* file_version */){
    05309             serialization::collection_size_type s1 (size1_);
    05310             serialization::collection_size_type s2 (size2_);
    05311             ar & serialization::make_nvp("size1",s1);
    05312             ar & serialization::make_nvp("size2",s2);
    05313             if (Archive::is_loading::value) {
    05314                 size1_ = s1;
    05315                 size2_ = s2;
    05316             }
    05317             ar & serialization::make_nvp("capacity", capacity_);
    05318             ar & serialization::make_nvp("filled", filled_);
    05319             ar & serialization::make_nvp("sorted_filled", sorted_filled_);
    05320             ar & serialization::make_nvp("sorted", sorted_);
    05321             ar & serialization::make_nvp("index1_data", index1_data_);
    05322             ar & serialization::make_nvp("index2_data", index2_data_);
    05323             ar & serialization::make_nvp("value_data", value_data_);
    05324             storage_invariants();
    05325         }
    05326
    05327     private:
    05328         void storage_invariants () const
    05329         {
    05330             BOOST_UBLAS_CHECK (capacity_ == index1_data_.size (), internal_logic ());
    05331             BOOST_UBLAS_CHECK (capacity_ == index2_data_.size (), internal_logic ());
    05332             BOOST_UBLAS_CHECK (capacity_ == value_data_.size (), internal_logic ());
    05333             BOOST_UBLAS_CHECK (filled_ <= capacity_, internal_logic ());
    05334             BOOST_UBLAS_CHECK (sorted_filled_ <= filled_, internal_logic ());
    05335             BOOST_UBLAS_CHECK (sorted_ == (sorted_filled_ == filled_), internal_logic ());
    05336         }
    05337
    05338         size_type size1_;
    05339         size_type size2_;
    05340         array_size_type capacity_;
    05341         mutable array_size_type filled_;
    05342         mutable array_size_type sorted_filled_;
    05343         mutable bool sorted_;
    05344         mutable index_array_type index1_data_;
    05345         mutable index_array_type index2_data_;
    05346         mutable value_array_type value_data_;
    05347         static const value_type zero_;
    05348
    05349         BOOST_UBLAS_INLINE
    05350         static size_type zero_based (size_type k_based_index) {
    05351             return k_based_index - IB;
    05352         }
    05353         BOOST_UBLAS_INLINE
    05354         static size_type k_based (size_type zero_based_index) {
    05355             return zero_based_index + IB;
    05356         }
    05357
    05358         friend class iterator1;
    05359         friend class iterator2;
    05360         friend class const_iterator1;
    05361         friend class const_iterator2;
    05362     };
    05363
    05364     template<class T, class L, std::size_t IB, class IA, class TA>
    05365     const typename coordinate_matrix<T, L, IB, IA, TA>::value_type coordinate_matrix<T, L, IB, IA, TA>::zero_ = value_type/*zero*/();
    05366
    05367 }}}
    05368
    05369 #endif