Boost.uBlas 1.49
Linear Algebra in C++: matrices, vectors and numeric algorithms

matrix_sparse.hpp

Go to the documentation of this file.
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