Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:43:19

0001 //
0002 //  Copyright (c) 2000-2002
0003 //  Joerg Walter, Mathias Koch
0004 //
0005 //  Distributed under the Boost Software License, Version 1.0. (See
0006 //  accompanying file LICENSE_1_0.txt or copy at
0007 //  http://www.boost.org/LICENSE_1_0.txt)
0008 //
0009 //  The authors gratefully acknowledge the support of
0010 //  GeNeSys mbH & Co. KG in producing this work.
0011 //
0012 
0013 #ifndef _BOOST_UBLAS_VECTOR_SPARSE_
0014 #define _BOOST_UBLAS_VECTOR_SPARSE_
0015 
0016 #include <boost/config.hpp>
0017 
0018 // In debug mode, MSCV enables iterator debugging, which additional checks are
0019 // executed for consistency. So, when two iterators are compared, it is tested
0020 // that they point to elements of the same container. If the check fails, then
0021 // the program is aborted.
0022 //
0023 // When matrices MVOV are traversed by column and then by row, the previous
0024 // check fails.
0025 //
0026 // MVOV::iterator2 iter2 = mvov.begin2();
0027 // for (; iter2 != mvov.end() ; iter2++) {
0028 //    MVOV::iterator1 iter1 = iter2.begin();
0029 //    .....
0030 // }
0031 //
0032 // These additional checks in iterators are disabled in this file, but their
0033 // status are restored at the end of file.
0034 // https://msdn.microsoft.com/en-us/library/hh697468.aspx
0035 #ifdef BOOST_MSVC
0036 #define _BACKUP_ITERATOR_DEBUG_LEVEL _ITERATOR_DEBUG_LEVEL
0037 #undef _ITERATOR_DEBUG_LEVEL
0038 #define _ITERATOR_DEBUG_LEVEL 0
0039 #endif
0040 
0041 #include <boost/numeric/ublas/storage_sparse.hpp>
0042 #include <boost/numeric/ublas/vector_expression.hpp>
0043 #include <boost/numeric/ublas/detail/vector_assign.hpp>
0044 #if BOOST_UBLAS_TYPE_CHECK
0045 #include <boost/numeric/ublas/vector.hpp>
0046 #endif
0047 
0048 // Iterators based on ideas of Jeremy Siek
0049 
0050 namespace boost { namespace numeric { namespace ublas {
0051 
0052 #ifdef BOOST_UBLAS_STRICT_VECTOR_SPARSE
0053 
0054     template<class V>
0055     class sparse_vector_element:
0056        public container_reference<V> {
0057     public:
0058         typedef V vector_type;
0059         typedef typename V::size_type size_type;
0060         typedef typename V::value_type value_type;
0061         typedef const value_type &const_reference;
0062         typedef value_type *pointer;
0063 
0064     private:
0065         // Proxied element operations
0066         void get_d () const {
0067             pointer p = (*this) ().find_element (i_);
0068             if (p)
0069                 d_ = *p;
0070             else
0071                 d_ = value_type/*zero*/();
0072         }
0073 
0074         void set (const value_type &s) const {
0075             pointer p = (*this) ().find_element (i_);
0076             if (!p)
0077                 (*this) ().insert_element (i_, s);
0078             else
0079                 *p = s;
0080         }
0081 
0082     public:
0083         // Construction and destruction
0084         sparse_vector_element (vector_type &v, size_type i):
0085             container_reference<vector_type> (v), i_ (i) {
0086         }
0087         BOOST_UBLAS_INLINE
0088         sparse_vector_element (const sparse_vector_element &p):
0089             container_reference<vector_type> (p), i_ (p.i_) {}
0090         BOOST_UBLAS_INLINE
0091         ~sparse_vector_element () {
0092         }
0093 
0094         // Assignment
0095         BOOST_UBLAS_INLINE
0096         sparse_vector_element &operator = (const sparse_vector_element &p) {
0097             // Overide the implict copy assignment
0098             p.get_d ();
0099             set (p.d_);
0100             return *this;
0101         }
0102         template<class D>
0103         BOOST_UBLAS_INLINE
0104         sparse_vector_element &operator = (const D &d) {
0105             set (d);
0106             return *this;
0107         }
0108         template<class D>
0109         BOOST_UBLAS_INLINE
0110         sparse_vector_element &operator += (const D &d) {
0111             get_d ();
0112             d_ += d;
0113             set (d_);
0114             return *this;
0115         }
0116         template<class D>
0117         BOOST_UBLAS_INLINE
0118         sparse_vector_element &operator -= (const D &d) {
0119             get_d ();
0120             d_ -= d;
0121             set (d_);
0122             return *this;
0123         }
0124         template<class D>
0125         BOOST_UBLAS_INLINE
0126         sparse_vector_element &operator *= (const D &d) {
0127             get_d ();
0128             d_ *= d;
0129             set (d_);
0130             return *this;
0131         }
0132         template<class D>
0133         BOOST_UBLAS_INLINE
0134         sparse_vector_element &operator /= (const D &d) {
0135             get_d ();
0136             d_ /= d;
0137             set (d_);
0138             return *this;
0139         }
0140 
0141         // Comparison
0142         template<class D>
0143         BOOST_UBLAS_INLINE
0144         bool operator == (const D &d) const {
0145             get_d ();
0146             return d_ == d;
0147         }
0148         template<class D>
0149         BOOST_UBLAS_INLINE
0150         bool operator != (const D &d) const {
0151             get_d ();
0152             return d_ != d;
0153         }
0154 
0155         // Conversion - weak link in proxy as d_ is not a perfect alias for the element
0156         BOOST_UBLAS_INLINE
0157         operator const_reference () const {
0158             get_d ();
0159             return d_;
0160         }
0161 
0162         // Conversion to reference - may be invalidated
0163         BOOST_UBLAS_INLINE
0164         value_type& ref () const {
0165             const pointer p = (*this) ().find_element (i_);
0166             if (!p)
0167                 return (*this) ().insert_element (i_, value_type/*zero*/());
0168             else
0169                 return *p;
0170         }
0171 
0172     private:
0173         size_type i_;
0174         mutable value_type d_;
0175     };
0176 
0177     /*
0178      * Generalise explicit reference access
0179      */
0180     namespace detail {
0181         template <class R>
0182         struct element_reference {
0183             typedef R& reference;
0184             static reference get_reference (reference r)
0185             {
0186                 return r;
0187             }
0188         };
0189         template <class V>
0190         struct element_reference<sparse_vector_element<V> > {
0191             typedef typename V::value_type& reference;
0192             static reference get_reference (const sparse_vector_element<V>& sve)
0193             {
0194                 return sve.ref ();
0195             }
0196         };
0197     }
0198     template <class VER>
0199     typename detail::element_reference<VER>::reference ref (VER& ver) {
0200         return detail::element_reference<VER>::get_reference (ver);
0201     }
0202     template <class VER>
0203     typename detail::element_reference<VER>::reference ref (const VER& ver) {
0204         return detail::element_reference<VER>::get_reference (ver);
0205     }
0206 
0207 
0208     template<class V>
0209     struct type_traits<sparse_vector_element<V> > {
0210         typedef typename V::value_type element_type;
0211         typedef type_traits<sparse_vector_element<V> > self_type;
0212         typedef typename type_traits<element_type>::value_type value_type;
0213         typedef typename type_traits<element_type>::const_reference const_reference;
0214         typedef sparse_vector_element<V> reference;
0215         typedef typename type_traits<element_type>::real_type real_type;
0216         typedef typename type_traits<element_type>::precision_type precision_type;
0217 
0218         static const unsigned plus_complexity = type_traits<element_type>::plus_complexity;
0219         static const unsigned multiplies_complexity = type_traits<element_type>::multiplies_complexity;
0220 
0221         static
0222         BOOST_UBLAS_INLINE
0223         real_type real (const_reference t) {
0224             return type_traits<element_type>::real (t);
0225         }
0226         static
0227         BOOST_UBLAS_INLINE
0228         real_type imag (const_reference t) {
0229             return type_traits<element_type>::imag (t);
0230         }
0231         static
0232         BOOST_UBLAS_INLINE
0233         value_type conj (const_reference t) {
0234             return type_traits<element_type>::conj (t);
0235         }
0236 
0237         static
0238         BOOST_UBLAS_INLINE
0239         real_type type_abs (const_reference t) {
0240             return type_traits<element_type>::type_abs (t);
0241         }
0242         static
0243         BOOST_UBLAS_INLINE
0244         value_type type_sqrt (const_reference t) {
0245             return type_traits<element_type>::type_sqrt (t);
0246         }
0247 
0248         static
0249         BOOST_UBLAS_INLINE
0250         real_type norm_1 (const_reference t) {
0251             return type_traits<element_type>::norm_1 (t);
0252         }
0253         static
0254         BOOST_UBLAS_INLINE
0255         real_type norm_2 (const_reference t) {
0256             return type_traits<element_type>::norm_2 (t);
0257         }
0258         static
0259         BOOST_UBLAS_INLINE
0260         real_type norm_inf (const_reference t) {
0261             return type_traits<element_type>::norm_inf (t);
0262         }
0263 
0264         static
0265         BOOST_UBLAS_INLINE
0266         bool equals (const_reference t1, const_reference t2) {
0267             return type_traits<element_type>::equals (t1, t2);
0268         }
0269     };
0270 
0271     template<class V1, class T2>
0272     struct promote_traits<sparse_vector_element<V1>, T2> {
0273         typedef typename promote_traits<typename sparse_vector_element<V1>::value_type, T2>::promote_type promote_type;
0274     };
0275     template<class T1, class V2>
0276     struct promote_traits<T1, sparse_vector_element<V2> > {
0277         typedef typename promote_traits<T1, typename sparse_vector_element<V2>::value_type>::promote_type promote_type;
0278     };
0279     template<class V1, class V2>
0280     struct promote_traits<sparse_vector_element<V1>, sparse_vector_element<V2> > {
0281         typedef typename promote_traits<typename sparse_vector_element<V1>::value_type,
0282                                         typename sparse_vector_element<V2>::value_type>::promote_type promote_type;
0283     };
0284 
0285 #endif
0286 
0287 
0288     /** \brief Index map based sparse vector
0289      *
0290      * A sparse vector of values of type T of variable size. The sparse storage type A can be 
0291      * \c std::map<size_t, T> or \c map_array<size_t, T>. This means that only non-zero elements
0292      * are effectively stored.
0293      *
0294      * For a \f$n\f$-dimensional sparse vector,  and 0 <= i < n the non-zero elements \f$v_i\f$ 
0295      * are mapped to consecutive elements of the associative container, i.e. for elements 
0296      * \f$k = v_{i_1}\f$ and \f$k + 1 = v_{i_2}\f$ of the container, holds \f$i_1 < i_2\f$.
0297      *
0298      * Supported parameters for the adapted array are \c map_array<std::size_t, T> and 
0299      * \c map_std<std::size_t, T>. The latter is equivalent to \c std::map<std::size_t, T>.
0300      *
0301      * \tparam T the type of object stored in the vector (like double, float, complex, etc...)
0302      * \tparam A the type of Storage array
0303      */
0304     template<class T, class A>
0305     class mapped_vector:
0306         public vector_container<mapped_vector<T, A> > {
0307 
0308         typedef T &true_reference;
0309         typedef T *pointer;
0310         typedef const T *const_pointer;
0311         typedef mapped_vector<T, A> self_type;
0312     public:
0313 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
0314         using vector_container<self_type>::operator ();
0315 #endif
0316         typedef typename A::size_type size_type;
0317         typedef typename A::difference_type difference_type;
0318         typedef T value_type;
0319         typedef A array_type;
0320         typedef const value_type &const_reference;
0321 #ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
0322         typedef typename detail::map_traits<A,T>::reference reference;
0323 #else
0324         typedef sparse_vector_element<self_type> reference;
0325 #endif
0326         typedef const vector_reference<const self_type> const_closure_type;
0327         typedef vector_reference<self_type> closure_type;
0328         typedef self_type vector_temporary_type;
0329         typedef sparse_tag storage_category;
0330 
0331         // Construction and destruction
0332         BOOST_UBLAS_INLINE
0333         mapped_vector ():
0334             vector_container<self_type> (),
0335             size_ (0), data_ () {}
0336         BOOST_UBLAS_INLINE
0337         mapped_vector (size_type size, size_type non_zeros = 0):
0338             vector_container<self_type> (),
0339             size_ (size), data_ () {
0340             detail::map_reserve (data(), restrict_capacity (non_zeros));
0341         }
0342         BOOST_UBLAS_INLINE
0343         mapped_vector (const mapped_vector &v):
0344             vector_container<self_type> (),
0345             size_ (v.size_), data_ (v.data_) {}
0346         template<class AE>
0347         BOOST_UBLAS_INLINE
0348         mapped_vector (const vector_expression<AE> &ae, size_type non_zeros = 0):
0349             vector_container<self_type> (),
0350             size_ (ae ().size ()), data_ () {
0351             detail::map_reserve (data(), restrict_capacity (non_zeros));
0352             vector_assign<scalar_assign> (*this, ae);
0353         }
0354 
0355         // Accessors
0356         BOOST_UBLAS_INLINE
0357         size_type size () const {
0358             return size_;
0359         }
0360         BOOST_UBLAS_INLINE
0361         size_type nnz_capacity () const {
0362             return detail::map_capacity (data ());
0363         }
0364         BOOST_UBLAS_INLINE
0365         size_type nnz () const {
0366             return data (). size ();
0367         }
0368 
0369         // Storage accessors
0370         BOOST_UBLAS_INLINE
0371         const array_type &data () const {
0372             return data_;
0373         }
0374         BOOST_UBLAS_INLINE
0375         array_type &data () {
0376             return data_;
0377         }
0378 
0379         // Resizing
0380     private:
0381         BOOST_UBLAS_INLINE
0382         size_type restrict_capacity (size_type non_zeros) const {
0383             non_zeros = (std::min) (non_zeros, size_);
0384             return non_zeros;
0385         }
0386     public:
0387         BOOST_UBLAS_INLINE
0388         void resize (size_type size, bool preserve = true) {
0389             size_ = size;
0390             if (preserve) {
0391                 data ().erase (data ().lower_bound(size_), data ().end());
0392             }
0393             else {
0394                 data ().clear ();
0395             }
0396         }
0397 
0398         // Reserving
0399         BOOST_UBLAS_INLINE
0400                 void reserve (size_type non_zeros, bool /*preserve*/ = true) {
0401             detail::map_reserve (data (), restrict_capacity (non_zeros));
0402         }
0403 
0404         // Element support
0405         BOOST_UBLAS_INLINE
0406         pointer find_element (size_type i) {
0407             return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i));
0408         }
0409         BOOST_UBLAS_INLINE
0410         const_pointer find_element (size_type i) const {
0411             const_subiterator_type it (data ().find (i));
0412             if (it == data ().end ())
0413                 return 0;
0414             BOOST_UBLAS_CHECK ((*it).first == i, internal_logic ());   // broken map
0415             return &(*it).second;
0416         }
0417 
0418         // Element access
0419         BOOST_UBLAS_INLINE
0420         const_reference operator () (size_type i) const {
0421             BOOST_UBLAS_CHECK (i < size_, bad_index ());
0422             const_subiterator_type it (data ().find (i));
0423             if (it == data ().end ())
0424                 return zero_;
0425             BOOST_UBLAS_CHECK ((*it).first == i, internal_logic ());   // broken map
0426             return (*it).second;
0427         }
0428         BOOST_UBLAS_INLINE
0429         true_reference ref (size_type i) {
0430             BOOST_UBLAS_CHECK (i < size_, bad_index ());
0431             std::pair<subiterator_type, bool> ii (data ().insert (typename array_type::value_type (i, value_type/*zero*/())));
0432             BOOST_UBLAS_CHECK ((ii.first)->first == i, internal_logic ());   // broken map
0433             return (ii.first)->second;
0434         }
0435         BOOST_UBLAS_INLINE
0436         reference operator () (size_type i) {
0437 #ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
0438             return ref (i);
0439 #else
0440             BOOST_UBLAS_CHECK (i < size_, bad_index ());
0441             return reference (*this, i);
0442 #endif
0443         }
0444 
0445         BOOST_UBLAS_INLINE
0446         const_reference operator [] (size_type i) const {
0447             return (*this) (i);
0448         }
0449         BOOST_UBLAS_INLINE
0450         reference operator [] (size_type i) {
0451             return (*this) (i);
0452         }
0453 
0454         // Element assignment
0455         BOOST_UBLAS_INLINE
0456         true_reference insert_element (size_type i, const_reference t) {
0457             std::pair<subiterator_type, bool> ii = data ().insert (typename array_type::value_type (i, t));
0458             BOOST_UBLAS_CHECK (ii.second, bad_index ());        // duplicate element
0459             BOOST_UBLAS_CHECK ((ii.first)->first == i, internal_logic ());   // broken map
0460             if (!ii.second)     // existing element
0461                 (ii.first)->second = t;
0462             return (ii.first)->second;
0463         }
0464         BOOST_UBLAS_INLINE
0465         void erase_element (size_type i) {
0466             subiterator_type it = data ().find (i);
0467             if (it == data ().end ())
0468                 return;
0469             data ().erase (it);
0470         }
0471 
0472         // Zeroing
0473         BOOST_UBLAS_INLINE
0474         void clear () {
0475             data ().clear ();
0476         }
0477 
0478         // Assignment
0479         BOOST_UBLAS_INLINE
0480         mapped_vector &operator = (const mapped_vector &v) {
0481             if (this != &v) {
0482                 size_ = v.size_;
0483                 data () = v.data ();
0484             }
0485             return *this;
0486         }
0487         template<class C>          // Container assignment without temporary
0488         BOOST_UBLAS_INLINE
0489         mapped_vector &operator = (const vector_container<C> &v) {
0490             resize (v ().size (), false);
0491             assign (v);
0492             return *this;
0493         }
0494         BOOST_UBLAS_INLINE
0495         mapped_vector &assign_temporary (mapped_vector &v) {
0496             swap (v);
0497             return *this;
0498         }
0499         template<class AE>
0500         BOOST_UBLAS_INLINE
0501         mapped_vector &operator = (const vector_expression<AE> &ae) {
0502             self_type temporary (ae, detail::map_capacity (data()));
0503             return assign_temporary (temporary);
0504         }
0505         template<class AE>
0506         BOOST_UBLAS_INLINE
0507         mapped_vector &assign (const vector_expression<AE> &ae) {
0508             vector_assign<scalar_assign> (*this, ae);
0509             return *this;
0510         }
0511 
0512         // Computed assignment
0513         template<class AE>
0514         BOOST_UBLAS_INLINE
0515         mapped_vector &operator += (const vector_expression<AE> &ae) {
0516             self_type temporary (*this + ae, detail::map_capacity (data()));
0517             return assign_temporary (temporary);
0518         }
0519         template<class C>          // Container assignment without temporary
0520         BOOST_UBLAS_INLINE
0521         mapped_vector &operator += (const vector_container<C> &v) {
0522             plus_assign (v);
0523             return *this;
0524         }
0525         template<class AE>
0526         BOOST_UBLAS_INLINE
0527         mapped_vector &plus_assign (const vector_expression<AE> &ae) {
0528             vector_assign<scalar_plus_assign> (*this, ae);
0529             return *this;
0530         }
0531         template<class AE>
0532         BOOST_UBLAS_INLINE
0533         mapped_vector &operator -= (const vector_expression<AE> &ae) {
0534             self_type temporary (*this - ae, detail::map_capacity (data()));
0535             return assign_temporary (temporary);
0536         }
0537         template<class C>          // Container assignment without temporary
0538         BOOST_UBLAS_INLINE
0539         mapped_vector &operator -= (const vector_container<C> &v) {
0540             minus_assign (v);
0541             return *this;
0542         }
0543         template<class AE>
0544         BOOST_UBLAS_INLINE
0545         mapped_vector &minus_assign (const vector_expression<AE> &ae) {
0546             vector_assign<scalar_minus_assign> (*this, ae);
0547             return *this;
0548         }
0549         template<class AT>
0550         BOOST_UBLAS_INLINE
0551         mapped_vector &operator *= (const AT &at) {
0552             vector_assign_scalar<scalar_multiplies_assign> (*this, at);
0553             return *this;
0554         }
0555         template<class AT>
0556         BOOST_UBLAS_INLINE
0557         mapped_vector &operator /= (const AT &at) {
0558             vector_assign_scalar<scalar_divides_assign> (*this, at);
0559             return *this;
0560         }
0561 
0562         // Swapping
0563         BOOST_UBLAS_INLINE
0564         void swap (mapped_vector &v) {
0565             if (this != &v) {
0566                 std::swap (size_, v.size_);
0567                 data ().swap (v.data ());
0568             }
0569         }
0570         BOOST_UBLAS_INLINE
0571         friend void swap (mapped_vector &v1, mapped_vector &v2) {
0572             v1.swap (v2);
0573         }
0574 
0575         // Iterator types
0576     private:
0577         // Use storage iterator
0578         typedef typename A::const_iterator const_subiterator_type;
0579         typedef typename A::iterator subiterator_type;
0580 
0581         BOOST_UBLAS_INLINE
0582         true_reference at_element (size_type i) {
0583             BOOST_UBLAS_CHECK (i < size_, bad_index ());
0584             subiterator_type it (data ().find (i));
0585             BOOST_UBLAS_CHECK (it != data ().end(), bad_index ());
0586             BOOST_UBLAS_CHECK ((*it).first == i, internal_logic ());   // broken map
0587             return it->second;
0588         }
0589 
0590     public:
0591         class const_iterator;
0592         class iterator;
0593 
0594         // Element lookup
0595         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
0596         const_iterator find (size_type i) const {
0597             return const_iterator (*this, data ().lower_bound (i));
0598         }
0599         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
0600         iterator find (size_type i) {
0601             return iterator (*this, data ().lower_bound (i));
0602         }
0603 
0604 
0605         class const_iterator:
0606             public container_const_reference<mapped_vector>,
0607             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
0608                                                const_iterator, value_type> {
0609         public:
0610             typedef typename mapped_vector::value_type value_type;
0611             typedef typename mapped_vector::difference_type difference_type;
0612             typedef typename mapped_vector::const_reference reference;
0613             typedef const typename mapped_vector::pointer pointer;
0614 
0615             // Construction and destruction
0616             BOOST_UBLAS_INLINE
0617             const_iterator ():
0618                 container_const_reference<self_type> (), it_ () {}
0619             BOOST_UBLAS_INLINE
0620             const_iterator (const self_type &v, const const_subiterator_type &it):
0621                 container_const_reference<self_type> (v), it_ (it) {}
0622             BOOST_UBLAS_INLINE
0623             const_iterator (const typename self_type::iterator &it):  // ISSUE self_type:: stops VC8 using std::iterator here
0624                 container_const_reference<self_type> (it ()), it_ (it.it_) {}
0625 
0626             // Arithmetic
0627             BOOST_UBLAS_INLINE
0628             const_iterator &operator ++ () {
0629                 ++ it_;
0630                 return *this;
0631             }
0632             BOOST_UBLAS_INLINE
0633             const_iterator &operator -- () {
0634                 -- it_;
0635                 return *this;
0636             }
0637 
0638             // Dereference
0639             BOOST_UBLAS_INLINE
0640             const_reference operator * () const {
0641                 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
0642                 return (*it_).second;
0643             }
0644 
0645             // Index
0646             BOOST_UBLAS_INLINE
0647             size_type index () const {
0648                 BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
0649                 BOOST_UBLAS_CHECK ((*it_).first < (*this) ().size (), bad_index ());
0650                 return (*it_).first;
0651             }
0652 
0653             // Assignment
0654             BOOST_UBLAS_INLINE
0655             const_iterator &operator = (const const_iterator &it) {
0656                 container_const_reference<self_type>::assign (&it ());
0657                 it_ = it.it_;
0658                 return *this;
0659             }
0660 
0661             // Comparison
0662             BOOST_UBLAS_INLINE
0663             bool operator == (const const_iterator &it) const {
0664                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
0665                 return it_ == it.it_;
0666             }
0667 
0668         private:
0669             const_subiterator_type it_;
0670         };
0671 
0672         BOOST_UBLAS_INLINE
0673         const_iterator begin () const {
0674             return const_iterator (*this, data ().begin ());
0675         }
0676         BOOST_UBLAS_INLINE
0677         const_iterator cbegin () const {
0678             return begin ();
0679         }
0680         BOOST_UBLAS_INLINE
0681         const_iterator end () const {
0682             return const_iterator (*this, data ().end ());
0683         }
0684         BOOST_UBLAS_INLINE
0685         const_iterator cend () const {
0686             return end ();
0687         }
0688 
0689         class iterator:
0690             public container_reference<mapped_vector>,
0691             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
0692                                                iterator, value_type> {
0693         public:
0694             typedef typename mapped_vector::value_type value_type;
0695             typedef typename mapped_vector::difference_type difference_type;
0696             typedef typename mapped_vector::true_reference reference;
0697             typedef typename mapped_vector::pointer pointer;
0698 
0699             // Construction and destruction
0700             BOOST_UBLAS_INLINE
0701             iterator ():
0702                 container_reference<self_type> (), it_ () {}
0703             BOOST_UBLAS_INLINE
0704             iterator (self_type &v, const subiterator_type &it):
0705                 container_reference<self_type> (v), it_ (it) {}
0706 
0707             // Arithmetic
0708             BOOST_UBLAS_INLINE
0709             iterator &operator ++ () {
0710                 ++ it_;
0711                 return *this;
0712             }
0713             BOOST_UBLAS_INLINE
0714             iterator &operator -- () {
0715                 -- it_;
0716                 return *this;
0717             }
0718 
0719             // Dereference
0720             BOOST_UBLAS_INLINE
0721             reference operator * () const {
0722                 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
0723                 return (*it_).second;
0724             }
0725 
0726             // Index
0727             BOOST_UBLAS_INLINE
0728             size_type index () const {
0729                 BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
0730                 BOOST_UBLAS_CHECK ((*it_).first < (*this) ().size (), bad_index ());
0731                 return (*it_).first;
0732             }
0733 
0734             // Assignment
0735             BOOST_UBLAS_INLINE
0736             iterator &operator = (const iterator &it) {
0737                 container_reference<self_type>::assign (&it ());
0738                 it_ = it.it_;
0739                 return *this;
0740             }
0741 
0742             // Comparison
0743             BOOST_UBLAS_INLINE
0744             bool operator == (const iterator &it) const {
0745                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
0746                 return it_ == it.it_;
0747             }
0748 
0749         private:
0750             subiterator_type it_;
0751 
0752             friend class const_iterator;
0753         };
0754 
0755         BOOST_UBLAS_INLINE
0756         iterator begin () {
0757             return iterator (*this, data ().begin ());
0758         }
0759         BOOST_UBLAS_INLINE
0760         iterator end () {
0761             return iterator (*this, data ().end ());
0762         }
0763 
0764         // Reverse iterator
0765         typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
0766         typedef reverse_iterator_base<iterator> reverse_iterator;
0767 
0768         BOOST_UBLAS_INLINE
0769         const_reverse_iterator rbegin () const {
0770             return const_reverse_iterator (end ());
0771         }
0772         BOOST_UBLAS_INLINE
0773         const_reverse_iterator crbegin () const {
0774             return rbegin ();
0775         }
0776         BOOST_UBLAS_INLINE
0777         const_reverse_iterator rend () const {
0778             return const_reverse_iterator (begin ());
0779         }
0780         BOOST_UBLAS_INLINE
0781         const_reverse_iterator crend () const {
0782             return rend ();
0783         }
0784         BOOST_UBLAS_INLINE
0785         reverse_iterator rbegin () {
0786             return reverse_iterator (end ());
0787         }
0788         BOOST_UBLAS_INLINE
0789         reverse_iterator rend () {
0790             return reverse_iterator (begin ());
0791         }
0792 
0793          // Serialization
0794         template<class Archive>
0795         void serialize(Archive & ar, const unsigned int /* file_version */){
0796             serialization::collection_size_type s (size_);
0797             ar & serialization::make_nvp("size",s);
0798             if (Archive::is_loading::value) {
0799                 size_ = s;
0800             }
0801             ar & serialization::make_nvp("data", data_);
0802         }
0803 
0804     private:
0805         size_type size_;
0806         array_type data_;
0807         static const value_type zero_;
0808     };
0809 
0810     template<class T, class A>
0811     const typename mapped_vector<T, A>::value_type mapped_vector<T, A>::zero_ = value_type/*zero*/();
0812 
0813 
0814     // Thanks to Kresimir Fresl for extending this to cover different index bases.
0815     
0816     /** \brief Compressed array based sparse vector
0817      *
0818      * a sparse vector of values of type T of variable size. The non zero values are stored as 
0819      * two seperate arrays: an index array and a value array. The index array is always sorted 
0820      * and there is at most one entry for each index. Inserting an element can be time consuming.
0821      * If the vector contains a few zero entries, then it is better to have a normal vector.
0822      * If the vector has a very high dimension with a few non-zero values, then this vector is
0823      * very memory efficient (at the cost of a few more computations).
0824      *
0825      * For a \f$n\f$-dimensional compressed vector and \f$0 \leq i < n\f$ the non-zero elements 
0826      * \f$v_i\f$ are mapped to consecutive elements of the index and value container, i.e. for 
0827      * elements \f$k = v_{i_1}\f$ and \f$k + 1 = v_{i_2}\f$ of these containers holds \f$i_1 < i_2\f$.
0828      *
0829      * Supported parameters for the adapted array (indices and values) are \c unbounded_array<> ,
0830      * \c bounded_array<> and \c std::vector<>.
0831      *
0832      * \tparam T the type of object stored in the vector (like double, float, complex, etc...)
0833      * \tparam IB the index base of the compressed vector. Default is 0. Other supported value is 1
0834      * \tparam IA the type of adapted array for indices. Default is \c unbounded_array<std::size_t>
0835      * \tparam TA the type of adapted array for values. Default is unbounded_array<T>
0836      */
0837     template<class T, std::size_t IB, class IA, class TA>
0838     class compressed_vector:
0839         public vector_container<compressed_vector<T, IB, IA, TA> > {
0840 
0841         typedef T &true_reference;
0842         typedef T *pointer;
0843         typedef const T *const_pointer;
0844         typedef compressed_vector<T, IB, IA, TA> self_type;
0845     public:
0846 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
0847         using vector_container<self_type>::operator ();
0848 #endif
0849         // ISSUE require type consistency check
0850         // is_convertable (IA::size_type, TA::size_type)
0851         typedef typename IA::value_type size_type;
0852         typedef typename IA::difference_type difference_type;
0853         typedef T value_type;
0854         typedef const T &const_reference;
0855 #ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
0856         typedef T &reference;
0857 #else
0858         typedef sparse_vector_element<self_type> reference;
0859 #endif
0860         typedef IA index_array_type;
0861         typedef TA value_array_type;
0862         typedef const vector_reference<const self_type> const_closure_type;
0863         typedef vector_reference<self_type> closure_type;
0864         typedef self_type vector_temporary_type;
0865         typedef sparse_tag storage_category;
0866 
0867         // Construction and destruction
0868         BOOST_UBLAS_INLINE
0869         compressed_vector ():
0870             vector_container<self_type> (),
0871             size_ (0), capacity_ (restrict_capacity (0)), filled_ (0),
0872             index_data_ (capacity_), value_data_ (capacity_) {
0873             storage_invariants ();
0874         }
0875         explicit BOOST_UBLAS_INLINE
0876         compressed_vector (size_type size, size_type non_zeros = 0):
0877             vector_container<self_type> (),
0878             size_ (size), capacity_ (restrict_capacity (non_zeros)), filled_ (0),
0879             index_data_ (capacity_), value_data_ (capacity_) {
0880         storage_invariants ();
0881         }
0882         BOOST_UBLAS_INLINE
0883         compressed_vector (const compressed_vector &v):
0884             vector_container<self_type> (),
0885             size_ (v.size_), capacity_ (v.capacity_), filled_ (v.filled_),
0886             index_data_ (v.index_data_), value_data_ (v.value_data_) {
0887             storage_invariants ();
0888         }
0889         template<class AE>
0890         BOOST_UBLAS_INLINE
0891         compressed_vector (const vector_expression<AE> &ae, size_type non_zeros = 0):
0892             vector_container<self_type> (),
0893             size_ (ae ().size ()), capacity_ (restrict_capacity (non_zeros)), filled_ (0),
0894             index_data_ (capacity_), value_data_ (capacity_) {
0895             storage_invariants ();
0896             vector_assign<scalar_assign> (*this, ae);
0897         }
0898 
0899         // Accessors
0900         BOOST_UBLAS_INLINE
0901         size_type size () const {
0902             return size_;
0903         }
0904         BOOST_UBLAS_INLINE
0905         size_type nnz_capacity () const {
0906             return capacity_;
0907         }
0908         BOOST_UBLAS_INLINE
0909         size_type nnz () const {
0910             return filled_;
0911         }
0912 
0913         // Storage accessors
0914         BOOST_UBLAS_INLINE
0915         static size_type index_base () {
0916             return IB;
0917         }
0918         BOOST_UBLAS_INLINE
0919         typename index_array_type::size_type filled () const {
0920             return filled_;
0921         }
0922         BOOST_UBLAS_INLINE
0923         const index_array_type &index_data () const {
0924             return index_data_;
0925         }
0926         BOOST_UBLAS_INLINE
0927         const value_array_type &value_data () const {
0928             return value_data_;
0929         }
0930         BOOST_UBLAS_INLINE
0931         void set_filled (const typename index_array_type::size_type & filled) {
0932             filled_ = filled;
0933             storage_invariants ();
0934         }
0935         BOOST_UBLAS_INLINE
0936         index_array_type &index_data () {
0937             return index_data_;
0938         }
0939         BOOST_UBLAS_INLINE
0940         value_array_type &value_data () {
0941             return value_data_;
0942         }
0943 
0944         // Resizing
0945     private:
0946         BOOST_UBLAS_INLINE
0947         size_type restrict_capacity (size_type non_zeros) const {
0948             non_zeros = (std::max) (non_zeros, size_type (1));
0949             non_zeros = (std::min) (non_zeros, size_);
0950             return non_zeros;
0951         }
0952     public:
0953         BOOST_UBLAS_INLINE
0954         void resize (size_type size, bool preserve = true) {
0955             size_ = size;
0956             capacity_ = restrict_capacity (capacity_);
0957             if (preserve) {
0958                 index_data_. resize (capacity_, size_type ());
0959                 value_data_. resize (capacity_, value_type ());
0960                 filled_ = (std::min) (capacity_, filled_);
0961                 while ((filled_ > 0) && (zero_based(index_data_[filled_ - 1]) >= size)) {
0962                     --filled_;
0963                 }
0964             }
0965             else {
0966                 index_data_. resize (capacity_);
0967                 value_data_. resize (capacity_);
0968                 filled_ = 0;
0969             }
0970             storage_invariants ();
0971         }
0972 
0973         // Reserving
0974         BOOST_UBLAS_INLINE
0975         void reserve (size_type non_zeros, bool preserve = true) {
0976             capacity_ = restrict_capacity (non_zeros);
0977             if (preserve) {
0978                 index_data_. resize (capacity_, size_type ());
0979                 value_data_. resize (capacity_, value_type ());
0980                 filled_ = (std::min) (capacity_, filled_);
0981             }
0982             else {
0983                 index_data_. resize (capacity_);
0984                 value_data_. resize (capacity_);
0985                 filled_ = 0;
0986             }
0987             storage_invariants ();
0988         }
0989 
0990         // Element support
0991         BOOST_UBLAS_INLINE
0992         pointer find_element (size_type i) {
0993             return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i));
0994         }
0995         BOOST_UBLAS_INLINE
0996         const_pointer find_element (size_type i) const {
0997             const_subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
0998             if (it == index_data_.begin () + filled_ || *it != k_based (i))
0999                 return 0;
1000             return &value_data_ [it - index_data_.begin ()];
1001         }
1002 
1003         // Element access
1004         BOOST_UBLAS_INLINE
1005         const_reference operator () (size_type i) const {
1006             BOOST_UBLAS_CHECK (i < size_, bad_index ());
1007             const_subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1008             if (it == index_data_.begin () + filled_ || *it != k_based (i))
1009                 return zero_;
1010             return value_data_ [it - index_data_.begin ()];
1011         }
1012         BOOST_UBLAS_INLINE
1013         true_reference ref (size_type i) {
1014             BOOST_UBLAS_CHECK (i < size_, bad_index ());
1015             subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1016             if (it == index_data_.begin () + filled_ || *it != k_based (i))
1017                 return insert_element (i, value_type/*zero*/());
1018             else
1019                 return value_data_ [it - index_data_.begin ()];
1020         }
1021         BOOST_UBLAS_INLINE
1022         reference operator () (size_type i) {
1023 #ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
1024             return ref (i) ;
1025 #else
1026             BOOST_UBLAS_CHECK (i < size_, bad_index ());
1027             return reference (*this, i);
1028 #endif
1029         }
1030 
1031         BOOST_UBLAS_INLINE
1032         const_reference operator [] (size_type i) const {
1033             return (*this) (i);
1034         }
1035         BOOST_UBLAS_INLINE
1036         reference operator [] (size_type i) {
1037             return (*this) (i);
1038         }
1039 
1040         // Element assignment
1041         BOOST_UBLAS_INLINE
1042         true_reference insert_element (size_type i, const_reference t) {
1043             BOOST_UBLAS_CHECK (!find_element (i), bad_index ());        // duplicate element
1044             if (filled_ >= capacity_)
1045                 reserve (2 * capacity_, true);
1046             subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1047             // ISSUE max_capacity limit due to difference_type
1048             typename std::iterator_traits<subiterator_type>::difference_type n = it - index_data_.begin ();
1049             BOOST_UBLAS_CHECK (filled_ == 0 || filled_ == typename index_array_type::size_type (n) || *it != k_based (i), internal_logic ());   // duplicate found by lower_bound
1050             ++ filled_;
1051             it = index_data_.begin () + n;
1052             std::copy_backward (it, index_data_.begin () + filled_ - 1, index_data_.begin () + filled_);
1053             *it = k_based (i);
1054             typename value_array_type::iterator itt (value_data_.begin () + n);
1055             std::copy_backward (itt, value_data_.begin () + filled_ - 1, value_data_.begin () + filled_);
1056             *itt = t;
1057             storage_invariants ();
1058             return *itt;
1059         }
1060         BOOST_UBLAS_INLINE
1061         void erase_element (size_type i) {
1062             subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1063             typename std::iterator_traits<subiterator_type>::difference_type  n = it - index_data_.begin ();
1064             if (filled_ > typename index_array_type::size_type (n) && *it == k_based (i)) {
1065                 std::copy (it + 1, index_data_.begin () + filled_, it);
1066                 typename value_array_type::iterator itt (value_data_.begin () + n);
1067                 std::copy (itt + 1, value_data_.begin () + filled_, itt);
1068                 -- filled_;
1069             }
1070             storage_invariants ();
1071         }
1072 
1073         // Zeroing
1074         BOOST_UBLAS_INLINE
1075         void clear () {
1076             filled_ = 0;
1077             storage_invariants ();
1078         }
1079 
1080         // Assignment
1081         BOOST_UBLAS_INLINE
1082         compressed_vector &operator = (const compressed_vector &v) {
1083             if (this != &v) {
1084                 size_ = v.size_;
1085                 capacity_ = v.capacity_;
1086                 filled_ = v.filled_;
1087                 index_data_ = v.index_data_;
1088                 value_data_ = v.value_data_;
1089             }
1090             storage_invariants ();
1091             return *this;
1092         }
1093         template<class C>          // Container assignment without temporary
1094         BOOST_UBLAS_INLINE
1095         compressed_vector &operator = (const vector_container<C> &v) {
1096             resize (v ().size (), false);
1097             assign (v);
1098             return *this;
1099         }
1100         BOOST_UBLAS_INLINE
1101         compressed_vector &assign_temporary (compressed_vector &v) {
1102             swap (v);
1103             return *this;
1104         }
1105         template<class AE>
1106         BOOST_UBLAS_INLINE
1107         compressed_vector &operator = (const vector_expression<AE> &ae) {
1108             self_type temporary (ae, capacity_);
1109             return assign_temporary (temporary);
1110         }
1111         template<class AE>
1112         BOOST_UBLAS_INLINE
1113         compressed_vector &assign (const vector_expression<AE> &ae) {
1114             vector_assign<scalar_assign> (*this, ae);
1115             return *this;
1116         }
1117 
1118         // Computed assignment
1119         template<class AE>
1120         BOOST_UBLAS_INLINE
1121         compressed_vector &operator += (const vector_expression<AE> &ae) {
1122             self_type temporary (*this + ae, capacity_);
1123             return assign_temporary (temporary);
1124         }
1125         template<class C>          // Container assignment without temporary
1126         BOOST_UBLAS_INLINE
1127         compressed_vector &operator += (const vector_container<C> &v) {
1128             plus_assign (v);
1129             return *this;
1130         }
1131         template<class AE>
1132         BOOST_UBLAS_INLINE
1133         compressed_vector &plus_assign (const vector_expression<AE> &ae) {
1134             vector_assign<scalar_plus_assign> (*this, ae);
1135             return *this;
1136         }
1137         template<class AE>
1138         BOOST_UBLAS_INLINE
1139         compressed_vector &operator -= (const vector_expression<AE> &ae) {
1140             self_type temporary (*this - ae, capacity_);
1141             return assign_temporary (temporary);
1142         }
1143         template<class C>          // Container assignment without temporary
1144         BOOST_UBLAS_INLINE
1145         compressed_vector &operator -= (const vector_container<C> &v) {
1146             minus_assign (v);
1147             return *this;
1148         }
1149         template<class AE>
1150         BOOST_UBLAS_INLINE
1151         compressed_vector &minus_assign (const vector_expression<AE> &ae) {
1152             vector_assign<scalar_minus_assign> (*this, ae);
1153             return *this;
1154         }
1155         template<class AT>
1156         BOOST_UBLAS_INLINE
1157         compressed_vector &operator *= (const AT &at) {
1158             vector_assign_scalar<scalar_multiplies_assign> (*this, at);
1159             return *this;
1160         }
1161         template<class AT>
1162         BOOST_UBLAS_INLINE
1163         compressed_vector &operator /= (const AT &at) {
1164             vector_assign_scalar<scalar_divides_assign> (*this, at);
1165             return *this;
1166         }
1167 
1168         // Swapping
1169         BOOST_UBLAS_INLINE
1170         void swap (compressed_vector &v) {
1171             if (this != &v) {
1172                 std::swap (size_, v.size_);
1173                 std::swap (capacity_, v.capacity_);
1174                 std::swap (filled_, v.filled_);
1175                 index_data_.swap (v.index_data_);
1176                 value_data_.swap (v.value_data_);
1177             }
1178             storage_invariants ();
1179         }
1180         BOOST_UBLAS_INLINE
1181         friend void swap (compressed_vector &v1, compressed_vector &v2) {
1182             v1.swap (v2);
1183         }
1184 
1185         // Back element insertion and erasure
1186         BOOST_UBLAS_INLINE
1187         void push_back (size_type i, const_reference t) {
1188             BOOST_UBLAS_CHECK (filled_ == 0 || index_data_ [filled_ - 1] < k_based (i), external_logic ());
1189             if (filled_ >= capacity_)
1190                 reserve (2 * capacity_, true);
1191             BOOST_UBLAS_CHECK (filled_ < capacity_, internal_logic ());
1192             index_data_ [filled_] = k_based (i);
1193             value_data_ [filled_] = t;
1194             ++ filled_;
1195             storage_invariants ();
1196         }
1197         BOOST_UBLAS_INLINE
1198         void pop_back () {
1199             BOOST_UBLAS_CHECK (filled_ > 0, external_logic ());
1200             -- filled_;
1201             storage_invariants ();
1202         }
1203 
1204         // Iterator types
1205     private:
1206         // Use index array iterator
1207         typedef typename IA::const_iterator const_subiterator_type;
1208         typedef typename IA::iterator subiterator_type;
1209 
1210         BOOST_UBLAS_INLINE
1211         true_reference at_element (size_type i) {
1212             BOOST_UBLAS_CHECK (i < size_, bad_index ());
1213             subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1214             BOOST_UBLAS_CHECK (it != index_data_.begin () + filled_ && *it == k_based (i), bad_index ());
1215             return value_data_ [it - index_data_.begin ()];
1216         }
1217 
1218     public:
1219         class const_iterator;
1220         class iterator;
1221 
1222         // Element lookup
1223         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1224         const_iterator find (size_type i) const {
1225             return const_iterator (*this, detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1226         }
1227         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1228         iterator find (size_type i) {
1229             return iterator (*this, detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1230         }
1231 
1232 
1233         class const_iterator:
1234             public container_const_reference<compressed_vector>,
1235             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
1236                                                const_iterator, value_type> {
1237         public:
1238             typedef typename compressed_vector::value_type value_type;
1239             typedef typename compressed_vector::difference_type difference_type;
1240             typedef typename compressed_vector::const_reference reference;
1241             typedef const typename compressed_vector::pointer pointer;
1242 
1243             // Construction and destruction
1244             BOOST_UBLAS_INLINE
1245             const_iterator ():
1246                 container_const_reference<self_type> (), it_ () {}
1247             BOOST_UBLAS_INLINE
1248             const_iterator (const self_type &v, const const_subiterator_type &it):
1249                 container_const_reference<self_type> (v), it_ (it) {}
1250             BOOST_UBLAS_INLINE
1251             const_iterator (const typename self_type::iterator &it):  // ISSUE self_type:: stops VC8 using std::iterator here
1252                 container_const_reference<self_type> (it ()), it_ (it.it_) {}
1253 
1254             // Arithmetic
1255             BOOST_UBLAS_INLINE
1256             const_iterator &operator ++ () {
1257                 ++ it_;
1258                 return *this;
1259             }
1260             BOOST_UBLAS_INLINE
1261             const_iterator &operator -- () {
1262                 -- it_;
1263                 return *this;
1264             }
1265 
1266             // Dereference
1267             BOOST_UBLAS_INLINE
1268             const_reference operator * () const {
1269                 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
1270                 return (*this) ().value_data_ [it_ - (*this) ().index_data_.begin ()];
1271             }
1272 
1273             // Index
1274             BOOST_UBLAS_INLINE
1275             size_type index () const {
1276                 BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
1277                 BOOST_UBLAS_CHECK ((*this) ().zero_based (*it_) < (*this) ().size (), bad_index ());
1278                 return (*this) ().zero_based (*it_);
1279             }
1280 
1281             // Assignment
1282             BOOST_UBLAS_INLINE
1283             const_iterator &operator = (const const_iterator &it) {
1284                 container_const_reference<self_type>::assign (&it ());
1285                 it_ = it.it_;
1286                 return *this;
1287             }
1288 
1289             // Comparison
1290             BOOST_UBLAS_INLINE
1291             bool operator == (const const_iterator &it) const {
1292                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1293                 return it_ == it.it_;
1294             }
1295 
1296         private:
1297             const_subiterator_type it_;
1298         };
1299 
1300         BOOST_UBLAS_INLINE
1301         const_iterator begin () const {
1302             return find (0);
1303         }
1304         BOOST_UBLAS_INLINE
1305         const_iterator cbegin () const {
1306             return begin ();
1307         }
1308         BOOST_UBLAS_INLINE
1309         const_iterator end () const {
1310             return find (size_);
1311         }
1312         BOOST_UBLAS_INLINE
1313         const_iterator cend () const {
1314             return end ();
1315         }
1316 
1317         class iterator:
1318             public container_reference<compressed_vector>,
1319             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
1320                                                iterator, value_type> {
1321         public:
1322             typedef typename compressed_vector::value_type value_type;
1323             typedef typename compressed_vector::difference_type difference_type;
1324             typedef typename compressed_vector::true_reference reference;
1325             typedef typename compressed_vector::pointer pointer;
1326 
1327             // Construction and destruction
1328             BOOST_UBLAS_INLINE
1329             iterator ():
1330                 container_reference<self_type> (), it_ () {}
1331             BOOST_UBLAS_INLINE
1332             iterator (self_type &v, const subiterator_type &it):
1333                 container_reference<self_type> (v), it_ (it) {}
1334 
1335             // Arithmetic
1336             BOOST_UBLAS_INLINE
1337             iterator &operator ++ () {
1338                 ++ it_;
1339                 return *this;
1340             }
1341             BOOST_UBLAS_INLINE
1342             iterator &operator -- () {
1343                 -- it_;
1344                 return *this;
1345             }
1346 
1347             // Dereference
1348             BOOST_UBLAS_INLINE
1349             reference operator * () const {
1350                 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
1351                 return (*this) ().value_data_ [it_ - (*this) ().index_data_.begin ()];
1352             }
1353 
1354             // Index
1355             BOOST_UBLAS_INLINE
1356             size_type index () const {
1357                 BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
1358                 BOOST_UBLAS_CHECK ((*this) ().zero_based (*it_) < (*this) ().size (), bad_index ());
1359                 return (*this) ().zero_based (*it_);
1360             }
1361 
1362             // Assignment
1363             BOOST_UBLAS_INLINE
1364             iterator &operator = (const iterator &it) {
1365                 container_reference<self_type>::assign (&it ());
1366                 it_ = it.it_;
1367                 return *this;
1368             }
1369 
1370             // Comparison
1371             BOOST_UBLAS_INLINE
1372             bool operator == (const iterator &it) const {
1373                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1374                 return it_ == it.it_;
1375             }
1376 
1377         private:
1378             subiterator_type it_;
1379 
1380             friend class const_iterator;
1381         };
1382 
1383         BOOST_UBLAS_INLINE
1384         iterator begin () {
1385             return find (0);
1386         }
1387         BOOST_UBLAS_INLINE
1388         iterator end () {
1389             return find (size_);
1390         }
1391 
1392         // Reverse iterator
1393         typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
1394         typedef reverse_iterator_base<iterator> reverse_iterator;
1395 
1396         BOOST_UBLAS_INLINE
1397         const_reverse_iterator rbegin () const {
1398             return const_reverse_iterator (end ());
1399         }
1400         BOOST_UBLAS_INLINE
1401         const_reverse_iterator crbegin () const {
1402             return rbegin ();
1403         }
1404         BOOST_UBLAS_INLINE
1405         const_reverse_iterator rend () const {
1406             return const_reverse_iterator (begin ());
1407         }
1408         BOOST_UBLAS_INLINE
1409         const_reverse_iterator crend () const {
1410             return rend ();
1411         }
1412         BOOST_UBLAS_INLINE
1413         reverse_iterator rbegin () {
1414             return reverse_iterator (end ());
1415         }
1416         BOOST_UBLAS_INLINE
1417         reverse_iterator rend () {
1418             return reverse_iterator (begin ());
1419         }
1420 
1421          // Serialization
1422         template<class Archive>
1423         void serialize(Archive & ar, const unsigned int /* file_version */){
1424             serialization::collection_size_type s (size_);
1425             ar & serialization::make_nvp("size",s);
1426             if (Archive::is_loading::value) {
1427                 size_ = s;
1428             }
1429             // ISSUE: filled may be much less than capacity
1430             // ISSUE: index_data_ and value_data_ are undefined between filled and capacity (trouble with 'nan'-values)
1431             ar & serialization::make_nvp("capacity", capacity_);
1432             ar & serialization::make_nvp("filled", filled_);
1433             ar & serialization::make_nvp("index_data", index_data_);
1434             ar & serialization::make_nvp("value_data", value_data_);
1435             storage_invariants();
1436         }
1437 
1438     private:
1439         void storage_invariants () const
1440         {
1441             BOOST_UBLAS_CHECK (capacity_ == index_data_.size (), internal_logic ());
1442             BOOST_UBLAS_CHECK (capacity_ == value_data_.size (), internal_logic ());
1443             BOOST_UBLAS_CHECK (filled_ <= capacity_, internal_logic ());
1444             BOOST_UBLAS_CHECK ((0 == filled_) || (zero_based(index_data_[filled_ - 1]) < size_), internal_logic ());
1445         }
1446 
1447         size_type size_;
1448         typename index_array_type::size_type capacity_;
1449         typename index_array_type::size_type filled_;
1450         index_array_type index_data_;
1451         value_array_type value_data_;
1452         static const value_type zero_;
1453 
1454         BOOST_UBLAS_INLINE
1455         static size_type zero_based (size_type k_based_index) {
1456             return k_based_index - IB;
1457         }
1458         BOOST_UBLAS_INLINE
1459         static size_type k_based (size_type zero_based_index) {
1460             return zero_based_index + IB;
1461         }
1462 
1463         friend class iterator;
1464         friend class const_iterator;
1465     };
1466 
1467     template<class T, std::size_t IB, class IA, class TA>
1468     const typename compressed_vector<T, IB, IA, TA>::value_type compressed_vector<T, IB, IA, TA>::zero_ = value_type/*zero*/();
1469 
1470     // Thanks to Kresimir Fresl for extending this to cover different index bases.
1471 
1472     /** \brief Coordimate array based sparse vector
1473      *
1474      * a sparse vector of values of type \c T of variable size. The non zero values are stored 
1475      * as two seperate arrays: an index array and a value array. The arrays may be out of order 
1476      * with multiple entries for each vector element. If there are multiple values for the same 
1477      * index the sum of these values is the real value. It is way more efficient for inserting values
1478      * than a \c compressed_vector but less memory efficient. Also linearly parsing a vector can 
1479      * be longer in specific cases than a \c compressed_vector.
1480      *
1481      * For a n-dimensional sorted coordinate vector and \f$ 0 \leq i < n\f$ the non-zero elements 
1482      * \f$v_i\f$ are mapped to consecutive elements of the index and value container, i.e. for 
1483      * elements \f$k = v_{i_1}\f$ and \f$k + 1 = v_{i_2}\f$ of these containers holds \f$i_1 < i_2\f$.
1484      *
1485      * Supported parameters for the adapted array (indices and values) are \c unbounded_array<> ,
1486      * \c bounded_array<> and \c std::vector<>.
1487      *
1488      * \tparam T the type of object stored in the vector (like double, float, complex, etc...)
1489      * \tparam IB the index base of the compressed vector. Default is 0. Other supported value is 1
1490      * \tparam IA the type of adapted array for indices. Default is \c unbounded_array<std::size_t>
1491      * \tparam TA the type of adapted array for values. Default is unbounded_array<T>
1492      */
1493     template<class T, std::size_t IB, class IA, class TA>
1494     class coordinate_vector:
1495         public vector_container<coordinate_vector<T, IB, IA, TA> > {
1496 
1497         typedef T &true_reference;
1498         typedef T *pointer;
1499         typedef const T *const_pointer;
1500         typedef coordinate_vector<T, IB, IA, TA> self_type;
1501     public:
1502 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
1503         using vector_container<self_type>::operator ();
1504 #endif
1505         // ISSUE require type consistency check
1506         // is_convertable (IA::size_type, TA::size_type)
1507         typedef typename IA::value_type size_type;
1508         typedef typename IA::difference_type difference_type;
1509         typedef T value_type;
1510         typedef const T &const_reference;
1511 #ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
1512         typedef T &reference;
1513 #else
1514         typedef sparse_vector_element<self_type> reference;
1515 #endif
1516         typedef IA index_array_type;
1517         typedef TA value_array_type;
1518         typedef const vector_reference<const self_type> const_closure_type;
1519         typedef vector_reference<self_type> closure_type;
1520         typedef self_type vector_temporary_type;
1521         typedef sparse_tag storage_category;
1522 
1523         // Construction and destruction
1524         BOOST_UBLAS_INLINE
1525         coordinate_vector ():
1526             vector_container<self_type> (),
1527             size_ (0), capacity_ (restrict_capacity (0)),
1528             filled_ (0), sorted_filled_ (filled_), sorted_ (true),
1529             index_data_ (capacity_), value_data_ (capacity_) {
1530             storage_invariants ();
1531         }
1532         explicit BOOST_UBLAS_INLINE
1533         coordinate_vector (size_type size, size_type non_zeros = 0):
1534             vector_container<self_type> (),
1535             size_ (size), capacity_ (restrict_capacity (non_zeros)),
1536             filled_ (0), sorted_filled_ (filled_), sorted_ (true),
1537             index_data_ (capacity_), value_data_ (capacity_) {
1538             storage_invariants ();
1539         }
1540         BOOST_UBLAS_INLINE
1541         coordinate_vector (const coordinate_vector &v):
1542             vector_container<self_type> (),
1543             size_ (v.size_), capacity_ (v.capacity_),
1544             filled_ (v.filled_), sorted_filled_ (v.sorted_filled_), sorted_ (v.sorted_),
1545             index_data_ (v.index_data_), value_data_ (v.value_data_) {
1546             storage_invariants ();
1547         }
1548         template<class AE>
1549         BOOST_UBLAS_INLINE
1550         coordinate_vector (const vector_expression<AE> &ae, size_type non_zeros = 0):
1551             vector_container<self_type> (),
1552             size_ (ae ().size ()), capacity_ (restrict_capacity (non_zeros)),
1553             filled_ (0), sorted_filled_ (filled_), sorted_ (true),
1554             index_data_ (capacity_), value_data_ (capacity_) {
1555             storage_invariants ();
1556             vector_assign<scalar_assign> (*this, ae);
1557         }
1558 
1559         // Accessors
1560         BOOST_UBLAS_INLINE
1561         size_type size () const {
1562             return size_;
1563         }
1564         BOOST_UBLAS_INLINE
1565         size_type nnz_capacity () const {
1566             return capacity_;
1567         }
1568         BOOST_UBLAS_INLINE
1569         size_type nnz () const {
1570             return filled_;
1571         }
1572 
1573         // Storage accessors
1574         BOOST_UBLAS_INLINE
1575         static size_type index_base () {
1576             return IB;
1577         }
1578         BOOST_UBLAS_INLINE
1579         typename index_array_type::size_type filled () const {
1580             return filled_;
1581         }
1582         BOOST_UBLAS_INLINE
1583         const index_array_type &index_data () const {
1584             return index_data_;
1585         }
1586         BOOST_UBLAS_INLINE
1587         const value_array_type &value_data () const {
1588             return value_data_;
1589         }
1590         BOOST_UBLAS_INLINE
1591         void set_filled (const typename index_array_type::size_type &sorted, const typename index_array_type::size_type &filled) {
1592             sorted_filled_ = sorted;
1593             filled_ = filled;
1594             storage_invariants ();
1595         }
1596         BOOST_UBLAS_INLINE
1597         index_array_type &index_data () {
1598             return index_data_;
1599         }
1600         BOOST_UBLAS_INLINE
1601         value_array_type &value_data () {
1602             return value_data_;
1603         }
1604 
1605         // Resizing
1606     private:
1607         BOOST_UBLAS_INLINE
1608         size_type restrict_capacity (size_type non_zeros) const {
1609              // minimum non_zeros
1610              non_zeros = (std::max) (non_zeros, size_type (1));
1611              // ISSUE no maximum as coordinate may contain inserted duplicates
1612              return non_zeros;
1613         }
1614     public:
1615         BOOST_UBLAS_INLINE
1616         void resize (size_type size, bool preserve = true) {
1617             if (preserve)
1618                 sort ();    // remove duplicate elements.
1619             size_ = size;
1620             capacity_ = restrict_capacity (capacity_);
1621             if (preserve) {
1622                 index_data_. resize (capacity_, size_type ());
1623                 value_data_. resize (capacity_, value_type ());
1624                 filled_ = (std::min) (capacity_, filled_);
1625                 while ((filled_ > 0) && (zero_based(index_data_[filled_ - 1]) >= size)) {
1626                     --filled_;
1627                 }
1628             }
1629             else {
1630                 index_data_. resize (capacity_);
1631                 value_data_. resize (capacity_);
1632                 filled_ = 0;
1633             }
1634             sorted_filled_ = filled_;
1635             storage_invariants ();
1636         }
1637         // Reserving
1638         BOOST_UBLAS_INLINE
1639         void reserve (size_type non_zeros, bool preserve = true) {
1640             if (preserve)
1641                 sort ();    // remove duplicate elements.
1642             capacity_ = restrict_capacity (non_zeros);
1643             if (preserve) {
1644                 index_data_. resize (capacity_, size_type ());
1645                 value_data_. resize (capacity_, value_type ());
1646                 filled_ = (std::min) (capacity_, filled_);
1647                 }
1648             else {
1649                 index_data_. resize (capacity_);
1650                 value_data_. resize (capacity_);
1651                 filled_ = 0;
1652             }
1653             sorted_filled_ = filled_;
1654             storage_invariants ();
1655         }
1656 
1657         // Element support
1658         BOOST_UBLAS_INLINE
1659         pointer find_element (size_type i) {
1660             return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i));
1661         }
1662         BOOST_UBLAS_INLINE
1663         const_pointer find_element (size_type i) const {
1664             sort ();
1665             const_subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1666             if (it == index_data_.begin () + filled_ || *it != k_based (i))
1667                 return 0;
1668             return &value_data_ [it - index_data_.begin ()];
1669         }
1670 
1671         // Element access
1672         BOOST_UBLAS_INLINE
1673         const_reference operator () (size_type i) const {
1674             BOOST_UBLAS_CHECK (i < size_, bad_index ());
1675             sort ();
1676             const_subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1677             if (it == index_data_.begin () + filled_ || *it != k_based (i))
1678                 return zero_;
1679             return value_data_ [it - index_data_.begin ()];
1680         }
1681         BOOST_UBLAS_INLINE
1682         true_reference ref (size_type i) {
1683             BOOST_UBLAS_CHECK (i < size_, bad_index ());
1684             sort ();
1685             subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1686             if (it == index_data_.begin () + filled_ || *it != k_based (i))
1687                 return insert_element (i, value_type/*zero*/());
1688             else
1689                 return value_data_ [it - index_data_.begin ()];
1690         }
1691         BOOST_UBLAS_INLINE
1692         reference operator () (size_type i) {
1693 #ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
1694             return ref (i);
1695 #else
1696             BOOST_UBLAS_CHECK (i < size_, bad_index ());
1697             return reference (*this, i);
1698 #endif
1699         }
1700 
1701         BOOST_UBLAS_INLINE
1702         const_reference operator [] (size_type i) const {
1703             return (*this) (i);
1704         }
1705         BOOST_UBLAS_INLINE
1706         reference operator [] (size_type i) {
1707             return (*this) (i);
1708         }
1709 
1710         // Element assignment
1711         BOOST_UBLAS_INLINE
1712         void append_element (size_type i, const_reference t) {
1713             if (filled_ >= capacity_)
1714                 reserve (2 * filled_, true);
1715             BOOST_UBLAS_CHECK (filled_ < capacity_, internal_logic ());
1716             index_data_ [filled_] = k_based (i);
1717             value_data_ [filled_] = t;
1718             ++ filled_;
1719             sorted_ = false;
1720             storage_invariants ();
1721         }
1722         BOOST_UBLAS_INLINE
1723         true_reference insert_element (size_type i, const_reference t) {
1724             BOOST_UBLAS_CHECK (!find_element (i), bad_index ());        // duplicate element
1725             append_element (i, t);
1726             return value_data_ [filled_ - 1];
1727         }
1728         BOOST_UBLAS_INLINE
1729         void erase_element (size_type i) {
1730             sort ();
1731             subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1732             typename std::iterator_traits<subiterator_type>::difference_type n = it - index_data_.begin ();
1733             if (filled_ > typename index_array_type::size_type (n) && *it == k_based (i)) {
1734                 std::copy (it + 1, index_data_.begin () + filled_, it);
1735                 typename value_array_type::iterator itt (value_data_.begin () + n);
1736                 std::copy (itt + 1, value_data_.begin () + filled_, itt);
1737                 -- filled_;
1738                 sorted_filled_ = filled_;
1739             }
1740             storage_invariants ();
1741         }
1742 
1743         // Zeroing
1744         BOOST_UBLAS_INLINE
1745         void clear () {
1746             filled_ = 0;
1747             sorted_filled_ = filled_;
1748             sorted_ = true;
1749             storage_invariants ();
1750         }
1751 
1752         // Assignment
1753         BOOST_UBLAS_INLINE
1754         coordinate_vector &operator = (const coordinate_vector &v) {
1755             if (this != &v) {
1756                 size_ = v.size_;
1757                 capacity_ = v.capacity_;
1758                 filled_ = v.filled_;
1759                 sorted_filled_ = v.sorted_filled_;
1760                 sorted_ = v.sorted_;
1761                 index_data_ = v.index_data_;
1762                 value_data_ = v.value_data_;
1763             }
1764             storage_invariants ();
1765             return *this;
1766         }
1767         template<class C>          // Container assignment without temporary
1768         BOOST_UBLAS_INLINE
1769         coordinate_vector &operator = (const vector_container<C> &v) {
1770             resize (v ().size (), false);
1771             assign (v);
1772             return *this;
1773         }
1774         BOOST_UBLAS_INLINE
1775         coordinate_vector &assign_temporary (coordinate_vector &v) {
1776             swap (v);
1777             return *this;
1778         }
1779         template<class AE>
1780         BOOST_UBLAS_INLINE
1781         coordinate_vector &operator = (const vector_expression<AE> &ae) {
1782             self_type temporary (ae, capacity_);
1783             return assign_temporary (temporary);
1784         }
1785         template<class AE>
1786         BOOST_UBLAS_INLINE
1787         coordinate_vector &assign (const vector_expression<AE> &ae) {
1788             vector_assign<scalar_assign> (*this, ae);
1789             return *this;
1790         }
1791 
1792         // Computed assignment
1793         template<class AE>
1794         BOOST_UBLAS_INLINE
1795         coordinate_vector &operator += (const vector_expression<AE> &ae) {
1796             self_type temporary (*this + ae, capacity_);
1797             return assign_temporary (temporary);
1798         }
1799         template<class C>          // Container assignment without temporary
1800         BOOST_UBLAS_INLINE
1801         coordinate_vector &operator += (const vector_container<C> &v) {
1802             plus_assign (v);
1803             return *this;
1804         }
1805         template<class AE>
1806         BOOST_UBLAS_INLINE
1807         coordinate_vector &plus_assign (const vector_expression<AE> &ae) {
1808             vector_assign<scalar_plus_assign> (*this, ae);
1809             return *this;
1810         }
1811         template<class AE>
1812         BOOST_UBLAS_INLINE
1813         coordinate_vector &operator -= (const vector_expression<AE> &ae) {
1814             self_type temporary (*this - ae, capacity_);
1815             return assign_temporary (temporary);
1816         }
1817         template<class C>          // Container assignment without temporary
1818         BOOST_UBLAS_INLINE
1819         coordinate_vector &operator -= (const vector_container<C> &v) {
1820             minus_assign (v);
1821             return *this;
1822         }
1823         template<class AE>
1824         BOOST_UBLAS_INLINE
1825         coordinate_vector &minus_assign (const vector_expression<AE> &ae) {
1826             vector_assign<scalar_minus_assign> (*this, ae);
1827             return *this;
1828         }
1829         template<class AT>
1830         BOOST_UBLAS_INLINE
1831         coordinate_vector &operator *= (const AT &at) {
1832             vector_assign_scalar<scalar_multiplies_assign> (*this, at);
1833             return *this;
1834         }
1835         template<class AT>
1836         BOOST_UBLAS_INLINE
1837         coordinate_vector &operator /= (const AT &at) {
1838             vector_assign_scalar<scalar_divides_assign> (*this, at);
1839             return *this;
1840         }
1841 
1842         // Swapping
1843         BOOST_UBLAS_INLINE
1844         void swap (coordinate_vector &v) {
1845             if (this != &v) {
1846                 std::swap (size_, v.size_);
1847                 std::swap (capacity_, v.capacity_);
1848                 std::swap (filled_, v.filled_);
1849                 std::swap (sorted_filled_, v.sorted_filled_);
1850                 std::swap (sorted_, v.sorted_);
1851                 index_data_.swap (v.index_data_);
1852                 value_data_.swap (v.value_data_);
1853             }
1854             storage_invariants ();
1855         }
1856         BOOST_UBLAS_INLINE
1857         friend void swap (coordinate_vector &v1, coordinate_vector &v2) {
1858             v1.swap (v2);
1859         }
1860 
1861         // replacement if STL lower bound algorithm for use of inplace_merge
1862         size_type lower_bound (size_type beg, size_type end, size_type target) const {
1863             while (end > beg) {
1864                 size_type mid = (beg + end) / 2;
1865                 if (index_data_[mid] < index_data_[target]) {
1866                     beg = mid + 1;
1867                 } else {
1868                     end = mid;
1869                 }
1870             }
1871             return beg;
1872         }
1873 
1874         // specialized replacement of STL inplace_merge to avoid compilation
1875         // problems with respect to the array_triple iterator
1876         void inplace_merge (size_type beg, size_type mid, size_type end) const {
1877             size_type len_lef = mid - beg;
1878             size_type len_rig = end - mid;
1879 
1880             if (len_lef == 1 && len_rig == 1) {
1881                 if (index_data_[mid] < index_data_[beg]) {
1882                     std::swap(index_data_[beg], index_data_[mid]);
1883                     std::swap(value_data_[beg], value_data_[mid]);
1884                 }
1885             } else if (len_lef > 0 && len_rig > 0) {
1886                 size_type lef_mid, rig_mid;
1887                 if (len_lef >= len_rig) {
1888                     lef_mid = (beg + mid) / 2;
1889                     rig_mid = lower_bound(mid, end, lef_mid);
1890                 } else {
1891                     rig_mid = (mid + end) / 2;
1892                     lef_mid = lower_bound(beg, mid, rig_mid);
1893                 }
1894                 std::rotate(&index_data_[0] + lef_mid, &index_data_[0] + mid, &index_data_[0] + rig_mid);
1895                 std::rotate(&value_data_[0] + lef_mid, &value_data_[0] + mid, &value_data_[0] + rig_mid);
1896 
1897                 size_type new_mid = lef_mid + rig_mid - mid;
1898                 inplace_merge(beg, lef_mid, new_mid);
1899                 inplace_merge(new_mid, rig_mid, end);
1900             }
1901         }
1902 
1903         // Sorting and summation of duplicates
1904         BOOST_UBLAS_INLINE
1905         void sort () const {
1906             if (! sorted_ && filled_ > 0) {
1907                 typedef index_pair_array<index_array_type, value_array_type> array_pair;
1908                 array_pair ipa (filled_, index_data_, value_data_);
1909 #ifndef BOOST_UBLAS_COO_ALWAYS_DO_FULL_SORT
1910                 const typename array_pair::iterator iunsorted = ipa.begin () + sorted_filled_;
1911                 // sort new elements and merge
1912                 std::sort (iunsorted, ipa.end ());
1913                 inplace_merge(0, sorted_filled_, filled_);
1914 #else
1915                 const typename array_pair::iterator iunsorted = ipa.begin ();
1916                 std::sort (iunsorted, ipa.end ());
1917 #endif
1918 
1919                 // sum duplicates with += and remove
1920                 size_type filled = 0;
1921                 for (size_type i = 1; i < filled_; ++ i) {
1922                     if (index_data_ [filled] != index_data_ [i]) {
1923                         ++ filled;
1924                         if (filled != i) {
1925                             index_data_ [filled] = index_data_ [i];
1926                             value_data_ [filled] = value_data_ [i];
1927                         }
1928                     } else {
1929                         value_data_ [filled] += value_data_ [i];
1930                     }
1931                 }
1932                 filled_ = filled + 1;
1933                 sorted_filled_ = filled_;
1934                 sorted_ = true;
1935                 storage_invariants ();
1936             }
1937         }
1938 
1939         // Back element insertion and erasure
1940         BOOST_UBLAS_INLINE
1941         void push_back (size_type i, const_reference t) {
1942             // must maintain sort order
1943             BOOST_UBLAS_CHECK (sorted_ && (filled_ == 0 || index_data_ [filled_ - 1] < k_based (i)), external_logic ());
1944             if (filled_ >= capacity_)
1945                 reserve (2 * filled_, true);
1946             BOOST_UBLAS_CHECK (filled_ < capacity_, internal_logic ());
1947             index_data_ [filled_] = k_based (i);
1948             value_data_ [filled_] = t;
1949             ++ filled_;
1950             sorted_filled_ = filled_;
1951             storage_invariants ();
1952         }
1953         BOOST_UBLAS_INLINE
1954         void pop_back () {
1955             // ISSUE invariants could be simpilfied if sorted required as precondition
1956             BOOST_UBLAS_CHECK (filled_ > 0, external_logic ());
1957             -- filled_;
1958             sorted_filled_ = (std::min) (sorted_filled_, filled_);
1959             sorted_ = sorted_filled_ = filled_;
1960             storage_invariants ();
1961         }
1962 
1963         // Iterator types
1964     private:
1965         // Use index array iterator
1966         typedef typename IA::const_iterator const_subiterator_type;
1967         typedef typename IA::iterator subiterator_type;
1968 
1969         BOOST_UBLAS_INLINE
1970         true_reference at_element (size_type i) {
1971             BOOST_UBLAS_CHECK (i < size_, bad_index ());
1972             sort ();
1973             subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1974             BOOST_UBLAS_CHECK (it != index_data_.begin () + filled_ && *it == k_based (i), bad_index ());
1975             return value_data_ [it - index_data_.begin ()];
1976         }
1977 
1978     public:
1979         class const_iterator;
1980         class iterator;
1981 
1982         // Element lookup
1983         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1984         const_iterator find (size_type i) const {
1985             sort ();
1986             return const_iterator (*this, detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1987         }
1988         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1989         iterator find (size_type i) {
1990             sort ();
1991             return iterator (*this, detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1992         }
1993 
1994 
1995         class const_iterator:
1996             public container_const_reference<coordinate_vector>,
1997             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
1998                                                const_iterator, value_type> {
1999         public:
2000             typedef typename coordinate_vector::value_type value_type;
2001             typedef typename coordinate_vector::difference_type difference_type;
2002             typedef typename coordinate_vector::const_reference reference;
2003             typedef const typename coordinate_vector::pointer pointer;
2004 
2005             // Construction and destruction
2006             BOOST_UBLAS_INLINE
2007             const_iterator ():
2008                 container_const_reference<self_type> (), it_ () {}
2009             BOOST_UBLAS_INLINE
2010             const_iterator (const self_type &v, const const_subiterator_type &it):
2011                 container_const_reference<self_type> (v), it_ (it) {}
2012             BOOST_UBLAS_INLINE
2013             const_iterator (const typename self_type::iterator &it):  // ISSUE self_type:: stops VC8 using std::iterator here
2014                 container_const_reference<self_type> (it ()), it_ (it.it_) {}
2015 
2016             // Arithmetic
2017             BOOST_UBLAS_INLINE
2018             const_iterator &operator ++ () {
2019                 ++ it_;
2020                 return *this;
2021             }
2022             BOOST_UBLAS_INLINE
2023             const_iterator &operator -- () {
2024                 -- it_;
2025                 return *this;
2026             }
2027 
2028             // Dereference
2029             BOOST_UBLAS_INLINE
2030             const_reference operator * () const {
2031                 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
2032                 return (*this) ().value_data_ [it_ - (*this) ().index_data_.begin ()];
2033             }
2034 
2035             // Index
2036             BOOST_UBLAS_INLINE
2037             size_type index () const {
2038                 BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
2039                 BOOST_UBLAS_CHECK ((*this) ().zero_based (*it_) < (*this) ().size (), bad_index ());
2040                 return (*this) ().zero_based (*it_);
2041             }
2042 
2043             // Assignment
2044             BOOST_UBLAS_INLINE
2045             const_iterator &operator = (const const_iterator &it) {
2046                 container_const_reference<self_type>::assign (&it ());
2047                 it_ = it.it_;
2048                 return *this;
2049             }
2050 
2051             // Comparison
2052             BOOST_UBLAS_INLINE
2053             bool operator == (const const_iterator &it) const {
2054                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2055                 return it_ == it.it_;
2056             }
2057 
2058         private:
2059             const_subiterator_type it_;
2060         };
2061 
2062         BOOST_UBLAS_INLINE
2063         const_iterator begin () const {
2064             return find (0);
2065         }
2066         BOOST_UBLAS_INLINE
2067         const_iterator cbegin () const {
2068             return begin();
2069         }
2070         BOOST_UBLAS_INLINE
2071         const_iterator end () const {
2072             return find (size_);
2073         }
2074         BOOST_UBLAS_INLINE
2075         const_iterator cend () const {
2076             return end();
2077         }
2078 
2079         class iterator:
2080             public container_reference<coordinate_vector>,
2081             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
2082                                                iterator, value_type> {
2083         public:
2084             typedef typename coordinate_vector::value_type value_type;
2085             typedef typename coordinate_vector::difference_type difference_type;
2086             typedef typename coordinate_vector::true_reference reference;
2087             typedef typename coordinate_vector::pointer pointer;
2088 
2089             // Construction and destruction
2090             BOOST_UBLAS_INLINE
2091             iterator ():
2092                 container_reference<self_type> (), it_ () {}
2093             BOOST_UBLAS_INLINE
2094             iterator (self_type &v, const subiterator_type &it):
2095                 container_reference<self_type> (v), it_ (it) {}
2096 
2097             // Arithmetic
2098             BOOST_UBLAS_INLINE
2099             iterator &operator ++ () {
2100                 ++ it_;
2101                 return *this;
2102             }
2103             BOOST_UBLAS_INLINE
2104             iterator &operator -- () {
2105                 -- it_;
2106                 return *this;
2107             }
2108 
2109             // Dereference
2110             BOOST_UBLAS_INLINE
2111             reference operator * () const {
2112                 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
2113                 return (*this) ().value_data_ [it_ - (*this) ().index_data_.begin ()];
2114             }
2115 
2116             // Index
2117             BOOST_UBLAS_INLINE
2118             size_type index () const {
2119                 BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
2120                 BOOST_UBLAS_CHECK ((*this) ().zero_based (*it_) < (*this) ().size (), bad_index ());
2121                 return (*this) ().zero_based (*it_);
2122             }
2123 
2124             // Assignment
2125             BOOST_UBLAS_INLINE
2126             iterator &operator = (const iterator &it) {
2127                 container_reference<self_type>::assign (&it ());
2128                 it_ = it.it_;
2129                 return *this;
2130             }
2131 
2132             // Comparison
2133             BOOST_UBLAS_INLINE
2134             bool operator == (const iterator &it) const {
2135                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2136                 return it_ == it.it_;
2137             }
2138 
2139         private:
2140             subiterator_type it_;
2141 
2142             friend class const_iterator;
2143         };
2144 
2145         BOOST_UBLAS_INLINE
2146         iterator begin () {
2147             return find (0);
2148         }
2149         BOOST_UBLAS_INLINE
2150         iterator end () {
2151             return find (size_);
2152         }
2153 
2154         // Reverse iterator
2155         typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
2156         typedef reverse_iterator_base<iterator> reverse_iterator;
2157 
2158         BOOST_UBLAS_INLINE
2159         const_reverse_iterator rbegin () const {
2160             return const_reverse_iterator (end ());
2161         }
2162         BOOST_UBLAS_INLINE
2163         const_reverse_iterator crbegin () const {
2164             return rbegin ();
2165         }
2166         BOOST_UBLAS_INLINE
2167         const_reverse_iterator rend () const {
2168             return const_reverse_iterator (begin ());
2169         }
2170         BOOST_UBLAS_INLINE
2171         const_reverse_iterator crend () const {
2172             return rend ();
2173         }
2174         BOOST_UBLAS_INLINE
2175         reverse_iterator rbegin () {
2176             return reverse_iterator (end ());
2177         }
2178         BOOST_UBLAS_INLINE
2179         reverse_iterator rend () {
2180             return reverse_iterator (begin ());
2181         }
2182 
2183          // Serialization
2184         template<class Archive>
2185         void serialize(Archive & ar, const unsigned int /* file_version */){
2186             serialization::collection_size_type s (size_);
2187             ar & serialization::make_nvp("size",s);
2188             if (Archive::is_loading::value) {
2189                 size_ = s;
2190             }
2191             // ISSUE: filled may be much less than capacity
2192             // ISSUE: index_data_ and value_data_ are undefined between filled and capacity (trouble with 'nan'-values)
2193             ar & serialization::make_nvp("capacity", capacity_);
2194             ar & serialization::make_nvp("filled", filled_);
2195             ar & serialization::make_nvp("sorted_filled", sorted_filled_);
2196             ar & serialization::make_nvp("sorted", sorted_);
2197             ar & serialization::make_nvp("index_data", index_data_);
2198             ar & serialization::make_nvp("value_data", value_data_);
2199             storage_invariants();
2200         }
2201 
2202     private:
2203         void storage_invariants () const
2204         {
2205             BOOST_UBLAS_CHECK (capacity_ == index_data_.size (), internal_logic ());
2206             BOOST_UBLAS_CHECK (capacity_ == value_data_.size (), internal_logic ());
2207             BOOST_UBLAS_CHECK (filled_ <= capacity_, internal_logic ());
2208             BOOST_UBLAS_CHECK (sorted_filled_ <= filled_, internal_logic ());
2209             BOOST_UBLAS_CHECK (sorted_ == (sorted_filled_ == filled_), internal_logic ());
2210             BOOST_UBLAS_CHECK ((0 == filled_) || (zero_based(index_data_[filled_ - 1]) < size_), internal_logic ());
2211         }
2212 
2213         size_type size_;
2214         size_type capacity_;
2215         mutable typename index_array_type::size_type filled_;
2216         mutable typename index_array_type::size_type sorted_filled_;
2217         mutable bool sorted_;
2218         mutable index_array_type index_data_;
2219         mutable value_array_type value_data_;
2220         static const value_type zero_;
2221 
2222         BOOST_UBLAS_INLINE
2223         static size_type zero_based (size_type k_based_index) {
2224             return k_based_index - IB;
2225         }
2226         BOOST_UBLAS_INLINE
2227         static size_type k_based (size_type zero_based_index) {
2228             return zero_based_index + IB;
2229         }
2230 
2231         friend class iterator;
2232         friend class const_iterator;
2233     };
2234 
2235     template<class T, std::size_t IB, class IA, class TA>
2236     const typename coordinate_vector<T, IB, IA, TA>::value_type coordinate_vector<T, IB, IA, TA>::zero_ = value_type/*zero*/();
2237 
2238 }}}
2239 
2240 #ifdef BOOST_MSVC
2241 #undef _ITERATOR_DEBUG_LEVEL
2242 #define _ITERATOR_DEBUG_LEVEL _BACKUP_ITERATOR_DEBUG_LEVEL
2243 #undef _BACKUP_ITERATOR_DEBUG_LEVEL
2244 #endif
2245 
2246 #endif