ceph/src/boost/libs/numeric/ublas/include/boost/numeric/ublas/vector_sparse.hpp

   1 //
   2 //  Copyright (c) 2000-2002
   3 //  Joerg Walter, Mathias Koch
   4 //
   5 //  Distributed under the Boost Software License, Version 1.0. (See
   6 //  accompanying file LICENSE_1_0.txt or copy at
   7 //  http://www.boost.org/LICENSE_1_0.txt)
   8 //
   9 //  The authors gratefully acknowledge the support of
  10 //  GeNeSys mbH & Co. KG in producing this work.
  11 //
  12
  13 #ifndef _BOOST_UBLAS_VECTOR_SPARSE_
  14 #define _BOOST_UBLAS_VECTOR_SPARSE_
  15
  16 #include <boost/numeric/ublas/storage_sparse.hpp>
  17 #include <boost/numeric/ublas/vector_expression.hpp>
  18 #include <boost/numeric/ublas/detail/vector_assign.hpp>
  19 #if BOOST_UBLAS_TYPE_CHECK
  20 #include <boost/numeric/ublas/vector.hpp>
  21 #endif
  22
  23 // Iterators based on ideas of Jeremy Siek
  24
  25 namespace boost { namespace numeric { namespace ublas {
  26
  27 #ifdef BOOST_UBLAS_STRICT_VECTOR_SPARSE
  28
  29     template<class V>
  30     class sparse_vector_element:
  31        public container_reference<V> {
  32     public:
  33         typedef V vector_type;
  34         typedef typename V::size_type size_type;
  35         typedef typename V::value_type value_type;
  36         typedef const value_type &const_reference;
  37         typedef value_type *pointer;
  38
  39     private:
  40         // Proxied element operations
  41         void get_d () const {
  42             pointer p = (*this) ().find_element (i_);
  43             if (p)
  44                 d_ = *p;
  45             else
  46                 d_ = value_type/*zero*/();
  47         }
  48
  49         void set (const value_type &s) const {
  50             pointer p = (*this) ().find_element (i_);
  51             if (!p)
  52                 (*this) ().insert_element (i_, s);
  53             else
  54                 *p = s;
  55         }
  56
  57     public:
  58         // Construction and destruction
  59         sparse_vector_element (vector_type &v, size_type i):
  60             container_reference<vector_type> (v), i_ (i) {
  61         }
  62         BOOST_UBLAS_INLINE
  63         sparse_vector_element (const sparse_vector_element &p):
  64             container_reference<vector_type> (p), i_ (p.i_) {}
  65         BOOST_UBLAS_INLINE
  66         ~sparse_vector_element () {
  67         }
  68
  69         // Assignment
  70         BOOST_UBLAS_INLINE
  71         sparse_vector_element &operator = (const sparse_vector_element &p) {
  72             // Overide the implict copy assignment
  73             p.get_d ();
  74             set (p.d_);
  75             return *this;
  76         }
  77         template<class D>
  78         BOOST_UBLAS_INLINE
  79         sparse_vector_element &operator = (const D &d) {
  80             set (d);
  81             return *this;
  82         }
  83         template<class D>
  84         BOOST_UBLAS_INLINE
  85         sparse_vector_element &operator += (const D &d) {
  86             get_d ();
  87             d_ += d;
  88             set (d_);
  89             return *this;
  90         }
  91         template<class D>
  92         BOOST_UBLAS_INLINE
  93         sparse_vector_element &operator -= (const D &d) {
  94             get_d ();
  95             d_ -= d;
  96             set (d_);
  97             return *this;
  98         }
  99         template<class D>
 100         BOOST_UBLAS_INLINE
 101         sparse_vector_element &operator *= (const D &d) {
 102             get_d ();
 103             d_ *= d;
 104             set (d_);
 105             return *this;
 106         }
 107         template<class D>
 108         BOOST_UBLAS_INLINE
 109         sparse_vector_element &operator /= (const D &d) {
 110             get_d ();
 111             d_ /= d;
 112             set (d_);
 113             return *this;
 114         }
 115
 116         // Comparison
 117         template<class D>
 118         BOOST_UBLAS_INLINE
 119         bool operator == (const D &d) const {
 120             get_d ();
 121             return d_ == d;
 122         }
 123         template<class D>
 124         BOOST_UBLAS_INLINE
 125         bool operator != (const D &d) const {
 126             get_d ();
 127             return d_ != d;
 128         }
 129
 130         // Conversion - weak link in proxy as d_ is not a perfect alias for the element
 131         BOOST_UBLAS_INLINE
 132         operator const_reference () const {
 133             get_d ();
 134             return d_;
 135         }
 136
 137         // Conversion to reference - may be invalidated
 138         BOOST_UBLAS_INLINE
 139         value_type& ref () const {
 140             const pointer p = (*this) ().find_element (i_);
 141             if (!p)
 142                 return (*this) ().insert_element (i_, value_type/*zero*/());
 143             else
 144                 return *p;
 145         }
 146
 147     private:
 148         size_type i_;
 149         mutable value_type d_;
 150     };
 151
 152     /*
 153      * Generalise explicit reference access
 154      */
 155     namespace detail {
 156         template <class R>
 157         struct element_reference {
 158             typedef R& reference;
 159             static reference get_reference (reference r)
 160             {
 161                 return r;
 162             }
 163         };
 164         template <class V>
 165         struct element_reference<sparse_vector_element<V> > {
 166             typedef typename V::value_type& reference;
 167             static reference get_reference (const sparse_vector_element<V>& sve)
 168             {
 169                 return sve.ref ();
 170             }
 171         };
 172     }
 173     template <class VER>
 174     typename detail::element_reference<VER>::reference ref (VER& ver) {
 175         return detail::element_reference<VER>::get_reference (ver);
 176     }
 177     template <class VER>
 178     typename detail::element_reference<VER>::reference ref (const VER& ver) {
 179         return detail::element_reference<VER>::get_reference (ver);
 180     }
 181
 182
 183     template<class V>
 184     struct type_traits<sparse_vector_element<V> > {
 185         typedef typename V::value_type element_type;
 186         typedef type_traits<sparse_vector_element<V> > self_type;
 187         typedef typename type_traits<element_type>::value_type value_type;
 188         typedef typename type_traits<element_type>::const_reference const_reference;
 189         typedef sparse_vector_element<V> reference;
 190         typedef typename type_traits<element_type>::real_type real_type;
 191         typedef typename type_traits<element_type>::precision_type precision_type;
 192
 193         static const unsigned plus_complexity = type_traits<element_type>::plus_complexity;
 194         static const unsigned multiplies_complexity = type_traits<element_type>::multiplies_complexity;
 195
 196         static
 197         BOOST_UBLAS_INLINE
 198         real_type real (const_reference t) {
 199             return type_traits<element_type>::real (t);
 200         }
 201         static
 202         BOOST_UBLAS_INLINE
 203         real_type imag (const_reference t) {
 204             return type_traits<element_type>::imag (t);
 205         }
 206         static
 207         BOOST_UBLAS_INLINE
 208         value_type conj (const_reference t) {
 209             return type_traits<element_type>::conj (t);
 210         }
 211
 212         static
 213         BOOST_UBLAS_INLINE
 214         real_type type_abs (const_reference t) {
 215             return type_traits<element_type>::type_abs (t);
 216         }
 217         static
 218         BOOST_UBLAS_INLINE
 219         value_type type_sqrt (const_reference t) {
 220             return type_traits<element_type>::type_sqrt (t);
 221         }
 222
 223         static
 224         BOOST_UBLAS_INLINE
 225         real_type norm_1 (const_reference t) {
 226             return type_traits<element_type>::norm_1 (t);
 227         }
 228         static
 229         BOOST_UBLAS_INLINE
 230         real_type norm_2 (const_reference t) {
 231             return type_traits<element_type>::norm_2 (t);
 232         }
 233         static
 234         BOOST_UBLAS_INLINE
 235         real_type norm_inf (const_reference t) {
 236             return type_traits<element_type>::norm_inf (t);
 237         }
 238
 239         static
 240         BOOST_UBLAS_INLINE
 241         bool equals (const_reference t1, const_reference t2) {
 242             return type_traits<element_type>::equals (t1, t2);
 243         }
 244     };
 245
 246     template<class V1, class T2>
 247     struct promote_traits<sparse_vector_element<V1>, T2> {
 248         typedef typename promote_traits<typename sparse_vector_element<V1>::value_type, T2>::promote_type promote_type;
 249     };
 250     template<class T1, class V2>
 251     struct promote_traits<T1, sparse_vector_element<V2> > {
 252         typedef typename promote_traits<T1, typename sparse_vector_element<V2>::value_type>::promote_type promote_type;
 253     };
 254     template<class V1, class V2>
 255     struct promote_traits<sparse_vector_element<V1>, sparse_vector_element<V2> > {
 256         typedef typename promote_traits<typename sparse_vector_element<V1>::value_type,
 257                                         typename sparse_vector_element<V2>::value_type>::promote_type promote_type;
 258     };
 259
 260 #endif
 261
 262
 263     /** \brief Index map based sparse vector
 264      *
 265      * A sparse vector of values of type T of variable size. The sparse storage type A can be
 266      * \c std::map<size_t, T> or \c map_array<size_t, T>. This means that only non-zero elements
 267      * are effectively stored.
 268      *
 269      * For a \f$n\f$-dimensional sparse vector,  and 0 <= i < n the non-zero elements \f$v_i\f$
 270      * are mapped to consecutive elements of the associative container, i.e. for elements
 271      * \f$k = v_{i_1}\f$ and \f$k + 1 = v_{i_2}\f$ of the container, holds \f$i_1 < i_2\f$.
 272      *
 273      * Supported parameters for the adapted array are \c map_array<std::size_t, T> and
 274      * \c map_std<std::size_t, T>. The latter is equivalent to \c std::map<std::size_t, T>.
 275      *
 276      * \tparam T the type of object stored in the vector (like double, float, complex, etc...)
 277      * \tparam A the type of Storage array
 278      */
 279     template<class T, class A>
 280     class mapped_vector:
 281         public vector_container<mapped_vector<T, A> > {
 282
 283         typedef T &true_reference;
 284         typedef T *pointer;
 285         typedef const T *const_pointer;
 286         typedef mapped_vector<T, A> self_type;
 287     public:
 288 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
 289         using vector_container<self_type>::operator ();
 290 #endif
 291         typedef typename A::size_type size_type;
 292         typedef typename A::difference_type difference_type;
 293         typedef T value_type;
 294         typedef A array_type;
 295         typedef const value_type &const_reference;
 296 #ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
 297         typedef typename detail::map_traits<A,T>::reference reference;
 298 #else
 299         typedef sparse_vector_element<self_type> reference;
 300 #endif
 301         typedef const vector_reference<const self_type> const_closure_type;
 302         typedef vector_reference<self_type> closure_type;
 303         typedef self_type vector_temporary_type;
 304         typedef sparse_tag storage_category;
 305
 306         // Construction and destruction
 307         BOOST_UBLAS_INLINE
 308         mapped_vector ():
 309             vector_container<self_type> (),
 310             size_ (0), data_ () {}
 311         BOOST_UBLAS_INLINE
 312         mapped_vector (size_type size, size_type non_zeros = 0):
 313             vector_container<self_type> (),
 314             size_ (size), data_ () {
 315             detail::map_reserve (data(), restrict_capacity (non_zeros));
 316         }
 317         BOOST_UBLAS_INLINE
 318         mapped_vector (const mapped_vector &v):
 319             vector_container<self_type> (),
 320             size_ (v.size_), data_ (v.data_) {}
 321         template<class AE>
 322         BOOST_UBLAS_INLINE
 323         mapped_vector (const vector_expression<AE> &ae, size_type non_zeros = 0):
 324             vector_container<self_type> (),
 325             size_ (ae ().size ()), data_ () {
 326             detail::map_reserve (data(), restrict_capacity (non_zeros));
 327             vector_assign<scalar_assign> (*this, ae);
 328         }
 329
 330         // Accessors
 331         BOOST_UBLAS_INLINE
 332         size_type size () const {
 333             return size_;
 334         }
 335         BOOST_UBLAS_INLINE
 336         size_type nnz_capacity () const {
 337             return detail::map_capacity (data ());
 338         }
 339         BOOST_UBLAS_INLINE
 340         size_type nnz () const {
 341             return data (). size ();
 342         }
 343
 344         // Storage accessors
 345         BOOST_UBLAS_INLINE
 346         const array_type &data () const {
 347             return data_;
 348         }
 349         BOOST_UBLAS_INLINE
 350         array_type &data () {
 351             return data_;
 352         }
 353
 354         // Resizing
 355     private:
 356         BOOST_UBLAS_INLINE
 357         size_type restrict_capacity (size_type non_zeros) const {
 358             non_zeros = (std::min) (non_zeros, size_);
 359             return non_zeros;
 360         }
 361     public:
 362         BOOST_UBLAS_INLINE
 363         void resize (size_type size, bool preserve = true) {
 364             size_ = size;
 365             if (preserve) {
 366                 data ().erase (data ().lower_bound(size_), data ().end());
 367             }
 368             else {
 369                 data ().clear ();
 370             }
 371         }
 372
 373         // Reserving
 374         BOOST_UBLAS_INLINE
 375         void reserve (size_type non_zeros = 0, bool preserve = true) {
 376             detail::map_reserve (data (), restrict_capacity (non_zeros));
 377         }
 378
 379         // Element support
 380         BOOST_UBLAS_INLINE
 381         pointer find_element (size_type i) {
 382             return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i));
 383         }
 384         BOOST_UBLAS_INLINE
 385         const_pointer find_element (size_type i) const {
 386             const_subiterator_type it (data ().find (i));
 387             if (it == data ().end ())
 388                 return 0;
 389             BOOST_UBLAS_CHECK ((*it).first == i, internal_logic ());   // broken map
 390             return &(*it).second;
 391         }
 392
 393         // Element access
 394         BOOST_UBLAS_INLINE
 395         const_reference operator () (size_type i) const {
 396             BOOST_UBLAS_CHECK (i < size_, bad_index ());
 397             const_subiterator_type it (data ().find (i));
 398             if (it == data ().end ())
 399                 return zero_;
 400             BOOST_UBLAS_CHECK ((*it).first == i, internal_logic ());   // broken map
 401             return (*it).second;
 402         }
 403         BOOST_UBLAS_INLINE
 404         true_reference ref (size_type i) {
 405             BOOST_UBLAS_CHECK (i < size_, bad_index ());
 406             std::pair<subiterator_type, bool> ii (data ().insert (typename array_type::value_type (i, value_type/*zero*/())));
 407             BOOST_UBLAS_CHECK ((ii.first)->first == i, internal_logic ());   // broken map
 408             return (ii.first)->second;
 409         }
 410         BOOST_UBLAS_INLINE
 411         reference operator () (size_type i) {
 412 #ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
 413             return ref (i);
 414 #else
 415             BOOST_UBLAS_CHECK (i < size_, bad_index ());
 416             return reference (*this, i);
 417 #endif
 418         }
 419
 420         BOOST_UBLAS_INLINE
 421         const_reference operator [] (size_type i) const {
 422             return (*this) (i);
 423         }
 424         BOOST_UBLAS_INLINE
 425         reference operator [] (size_type i) {
 426             return (*this) (i);
 427         }
 428
 429         // Element assignment
 430         BOOST_UBLAS_INLINE
 431         true_reference insert_element (size_type i, const_reference t) {
 432             std::pair<subiterator_type, bool> ii = data ().insert (typename array_type::value_type (i, t));
 433             BOOST_UBLAS_CHECK (ii.second, bad_index ());        // duplicate element
 434             BOOST_UBLAS_CHECK ((ii.first)->first == i, internal_logic ());   // broken map
 435             if (!ii.second)     // existing element
 436                 (ii.first)->second = t;
 437             return (ii.first)->second;
 438         }
 439         BOOST_UBLAS_INLINE
 440         void erase_element (size_type i) {
 441             subiterator_type it = data ().find (i);
 442             if (it == data ().end ())
 443                 return;
 444             data ().erase (it);
 445         }
 446
 447         // Zeroing
 448         BOOST_UBLAS_INLINE
 449         void clear () {
 450             data ().clear ();
 451         }
 452
 453         // Assignment
 454         BOOST_UBLAS_INLINE
 455         mapped_vector &operator = (const mapped_vector &v) {
 456             if (this != &v) {
 457                 size_ = v.size_;
 458                 data () = v.data ();
 459             }
 460             return *this;
 461         }
 462         template<class C>          // Container assignment without temporary
 463         BOOST_UBLAS_INLINE
 464         mapped_vector &operator = (const vector_container<C> &v) {
 465             resize (v ().size (), false);
 466             assign (v);
 467             return *this;
 468         }
 469         BOOST_UBLAS_INLINE
 470         mapped_vector &assign_temporary (mapped_vector &v) {
 471             swap (v);
 472             return *this;
 473         }
 474         template<class AE>
 475         BOOST_UBLAS_INLINE
 476         mapped_vector &operator = (const vector_expression<AE> &ae) {
 477             self_type temporary (ae, detail::map_capacity (data()));
 478             return assign_temporary (temporary);
 479         }
 480         template<class AE>
 481         BOOST_UBLAS_INLINE
 482         mapped_vector &assign (const vector_expression<AE> &ae) {
 483             vector_assign<scalar_assign> (*this, ae);
 484             return *this;
 485         }
 486
 487         // Computed assignment
 488         template<class AE>
 489         BOOST_UBLAS_INLINE
 490         mapped_vector &operator += (const vector_expression<AE> &ae) {
 491             self_type temporary (*this + ae, detail::map_capacity (data()));
 492             return assign_temporary (temporary);
 493         }
 494         template<class C>          // Container assignment without temporary
 495         BOOST_UBLAS_INLINE
 496         mapped_vector &operator += (const vector_container<C> &v) {
 497             plus_assign (v);
 498             return *this;
 499         }
 500         template<class AE>
 501         BOOST_UBLAS_INLINE
 502         mapped_vector &plus_assign (const vector_expression<AE> &ae) {
 503             vector_assign<scalar_plus_assign> (*this, ae);
 504             return *this;
 505         }
 506         template<class AE>
 507         BOOST_UBLAS_INLINE
 508         mapped_vector &operator -= (const vector_expression<AE> &ae) {
 509             self_type temporary (*this - ae, detail::map_capacity (data()));
 510             return assign_temporary (temporary);
 511         }
 512         template<class C>          // Container assignment without temporary
 513         BOOST_UBLAS_INLINE
 514         mapped_vector &operator -= (const vector_container<C> &v) {
 515             minus_assign (v);
 516             return *this;
 517         }
 518         template<class AE>
 519         BOOST_UBLAS_INLINE
 520         mapped_vector &minus_assign (const vector_expression<AE> &ae) {
 521             vector_assign<scalar_minus_assign> (*this, ae);
 522             return *this;
 523         }
 524         template<class AT>
 525         BOOST_UBLAS_INLINE
 526         mapped_vector &operator *= (const AT &at) {
 527             vector_assign_scalar<scalar_multiplies_assign> (*this, at);
 528             return *this;
 529         }
 530         template<class AT>
 531         BOOST_UBLAS_INLINE
 532         mapped_vector &operator /= (const AT &at) {
 533             vector_assign_scalar<scalar_divides_assign> (*this, at);
 534             return *this;
 535         }
 536
 537         // Swapping
 538         BOOST_UBLAS_INLINE
 539         void swap (mapped_vector &v) {
 540             if (this != &v) {
 541                 std::swap (size_, v.size_);
 542                 data ().swap (v.data ());
 543             }
 544         }
 545         BOOST_UBLAS_INLINE
 546         friend void swap (mapped_vector &v1, mapped_vector &v2) {
 547             v1.swap (v2);
 548         }
 549
 550         // Iterator types
 551     private:
 552         // Use storage iterator
 553         typedef typename A::const_iterator const_subiterator_type;
 554         typedef typename A::iterator subiterator_type;
 555
 556         BOOST_UBLAS_INLINE
 557         true_reference at_element (size_type i) {
 558             BOOST_UBLAS_CHECK (i < size_, bad_index ());
 559             subiterator_type it (data ().find (i));
 560             BOOST_UBLAS_CHECK (it != data ().end(), bad_index ());
 561             BOOST_UBLAS_CHECK ((*it).first == i, internal_logic ());   // broken map
 562             return it->second;
 563         }
 564
 565     public:
 566         class const_iterator;
 567         class iterator;
 568
 569         // Element lookup
 570         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
 571         const_iterator find (size_type i) const {
 572             return const_iterator (*this, data ().lower_bound (i));
 573         }
 574         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
 575         iterator find (size_type i) {
 576             return iterator (*this, data ().lower_bound (i));
 577         }
 578
 579
 580         class const_iterator:
 581             public container_const_reference<mapped_vector>,
 582             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
 583                                                const_iterator, value_type> {
 584         public:
 585             typedef typename mapped_vector::value_type value_type;
 586             typedef typename mapped_vector::difference_type difference_type;
 587             typedef typename mapped_vector::const_reference reference;
 588             typedef const typename mapped_vector::pointer pointer;
 589
 590             // Construction and destruction
 591             BOOST_UBLAS_INLINE
 592             const_iterator ():
 593                 container_const_reference<self_type> (), it_ () {}
 594             BOOST_UBLAS_INLINE
 595             const_iterator (const self_type &v, const const_subiterator_type &it):
 596                 container_const_reference<self_type> (v), it_ (it) {}
 597             BOOST_UBLAS_INLINE
 598             const_iterator (const typename self_type::iterator &it):  // ISSUE self_type:: stops VC8 using std::iterator here
 599                 container_const_reference<self_type> (it ()), it_ (it.it_) {}
 600
 601             // Arithmetic
 602             BOOST_UBLAS_INLINE
 603             const_iterator &operator ++ () {
 604                 ++ it_;
 605                 return *this;
 606             }
 607             BOOST_UBLAS_INLINE
 608             const_iterator &operator -- () {
 609                 -- it_;
 610                 return *this;
 611             }
 612
 613             // Dereference
 614             BOOST_UBLAS_INLINE
 615             const_reference operator * () const {
 616                 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
 617                 return (*it_).second;
 618             }
 619
 620             // Index
 621             BOOST_UBLAS_INLINE
 622             size_type index () const {
 623                 BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
 624                 BOOST_UBLAS_CHECK ((*it_).first < (*this) ().size (), bad_index ());
 625                 return (*it_).first;
 626             }
 627
 628             // Assignment
 629             BOOST_UBLAS_INLINE
 630             const_iterator &operator = (const const_iterator &it) {
 631                 container_const_reference<self_type>::assign (&it ());
 632                 it_ = it.it_;
 633                 return *this;
 634             }
 635
 636             // Comparison
 637             BOOST_UBLAS_INLINE
 638             bool operator == (const const_iterator &it) const {
 639                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
 640                 return it_ == it.it_;
 641             }
 642
 643         private:
 644             const_subiterator_type it_;
 645         };
 646
 647         BOOST_UBLAS_INLINE
 648         const_iterator begin () const {
 649             return const_iterator (*this, data ().begin ());
 650         }
 651         BOOST_UBLAS_INLINE
 652         const_iterator cbegin () const {
 653             return begin ();
 654         }
 655         BOOST_UBLAS_INLINE
 656         const_iterator end () const {
 657             return const_iterator (*this, data ().end ());
 658         }
 659         BOOST_UBLAS_INLINE
 660         const_iterator cend () const {
 661             return end ();
 662         }
 663
 664         class iterator:
 665             public container_reference<mapped_vector>,
 666             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
 667                                                iterator, value_type> {
 668         public:
 669             typedef typename mapped_vector::value_type value_type;
 670             typedef typename mapped_vector::difference_type difference_type;
 671             typedef typename mapped_vector::true_reference reference;
 672             typedef typename mapped_vector::pointer pointer;
 673
 674             // Construction and destruction
 675             BOOST_UBLAS_INLINE
 676             iterator ():
 677                 container_reference<self_type> (), it_ () {}
 678             BOOST_UBLAS_INLINE
 679             iterator (self_type &v, const subiterator_type &it):
 680                 container_reference<self_type> (v), it_ (it) {}
 681
 682             // Arithmetic
 683             BOOST_UBLAS_INLINE
 684             iterator &operator ++ () {
 685                 ++ it_;
 686                 return *this;
 687             }
 688             BOOST_UBLAS_INLINE
 689             iterator &operator -- () {
 690                 -- it_;
 691                 return *this;
 692             }
 693
 694             // Dereference
 695             BOOST_UBLAS_INLINE
 696             reference operator * () const {
 697                 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
 698                 return (*it_).second;
 699             }
 700
 701             // Index
 702             BOOST_UBLAS_INLINE
 703             size_type index () const {
 704                 BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
 705                 BOOST_UBLAS_CHECK ((*it_).first < (*this) ().size (), bad_index ());
 706                 return (*it_).first;
 707             }
 708
 709             // Assignment
 710             BOOST_UBLAS_INLINE
 711             iterator &operator = (const iterator &it) {
 712                 container_reference<self_type>::assign (&it ());
 713                 it_ = it.it_;
 714                 return *this;
 715             }
 716
 717             // Comparison
 718             BOOST_UBLAS_INLINE
 719             bool operator == (const iterator &it) const {
 720                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
 721                 return it_ == it.it_;
 722             }
 723
 724         private:
 725             subiterator_type it_;
 726
 727             friend class const_iterator;
 728         };
 729
 730         BOOST_UBLAS_INLINE
 731         iterator begin () {
 732             return iterator (*this, data ().begin ());
 733         }
 734         BOOST_UBLAS_INLINE
 735         iterator end () {
 736             return iterator (*this, data ().end ());
 737         }
 738
 739         // Reverse iterator
 740         typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
 741         typedef reverse_iterator_base<iterator> reverse_iterator;
 742
 743         BOOST_UBLAS_INLINE
 744         const_reverse_iterator rbegin () const {
 745             return const_reverse_iterator (end ());
 746         }
 747         BOOST_UBLAS_INLINE
 748         const_reverse_iterator crbegin () const {
 749             return rbegin ();
 750         }
 751         BOOST_UBLAS_INLINE
 752         const_reverse_iterator rend () const {
 753             return const_reverse_iterator (begin ());
 754         }
 755         BOOST_UBLAS_INLINE
 756         const_reverse_iterator crend () const {
 757             return rend ();
 758         }
 759         BOOST_UBLAS_INLINE
 760         reverse_iterator rbegin () {
 761             return reverse_iterator (end ());
 762         }
 763         BOOST_UBLAS_INLINE
 764         reverse_iterator rend () {
 765             return reverse_iterator (begin ());
 766         }
 767
 768          // Serialization
 769         template<class Archive>
 770         void serialize(Archive & ar, const unsigned int /* file_version */){
 771             serialization::collection_size_type s (size_);
 772             ar & serialization::make_nvp("size",s);
 773             if (Archive::is_loading::value) {
 774                 size_ = s;
 775             }
 776             ar & serialization::make_nvp("data", data_);
 777         }
 778
 779     private:
 780         size_type size_;
 781         array_type data_;
 782         static const value_type zero_;
 783     };
 784
 785     template<class T, class A>
 786     const typename mapped_vector<T, A>::value_type mapped_vector<T, A>::zero_ = value_type/*zero*/();
 787
 788
 789     // Thanks to Kresimir Fresl for extending this to cover different index bases.
 790
 791     /** \brief Compressed array based sparse vector
 792      *
 793      * a sparse vector of values of type T of variable size. The non zero values are stored as
 794      * two seperate arrays: an index array and a value array. The index array is always sorted
 795      * and there is at most one entry for each index. Inserting an element can be time consuming.
 796      * If the vector contains a few zero entries, then it is better to have a normal vector.
 797      * If the vector has a very high dimension with a few non-zero values, then this vector is
 798      * very memory efficient (at the cost of a few more computations).
 799      *
 800      * For a \f$n\f$-dimensional compressed vector and \f$0 \leq i < n\f$ the non-zero elements
 801      * \f$v_i\f$ are mapped to consecutive elements of the index and value container, i.e. for
 802      * 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$.
 803      *
 804      * Supported parameters for the adapted array (indices and values) are \c unbounded_array<> ,
 805      * \c bounded_array<> and \c std::vector<>.
 806      *
 807      * \tparam T the type of object stored in the vector (like double, float, complex, etc...)
 808      * \tparam IB the index base of the compressed vector. Default is 0. Other supported value is 1
 809      * \tparam IA the type of adapted array for indices. Default is \c unbounded_array<std::size_t>
 810      * \tparam TA the type of adapted array for values. Default is unbounded_array<T>
 811      */
 812     template<class T, std::size_t IB, class IA, class TA>
 813     class compressed_vector:
 814         public vector_container<compressed_vector<T, IB, IA, TA> > {
 815
 816         typedef T &true_reference;
 817         typedef T *pointer;
 818         typedef const T *const_pointer;
 819         typedef compressed_vector<T, IB, IA, TA> self_type;
 820     public:
 821 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
 822         using vector_container<self_type>::operator ();
 823 #endif
 824         // ISSUE require type consistency check
 825         // is_convertable (IA::size_type, TA::size_type)
 826         typedef typename IA::value_type size_type;
 827         typedef typename IA::difference_type difference_type;
 828         typedef T value_type;
 829         typedef const T &const_reference;
 830 #ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
 831         typedef T &reference;
 832 #else
 833         typedef sparse_vector_element<self_type> reference;
 834 #endif
 835         typedef IA index_array_type;
 836         typedef TA value_array_type;
 837         typedef const vector_reference<const self_type> const_closure_type;
 838         typedef vector_reference<self_type> closure_type;
 839         typedef self_type vector_temporary_type;
 840         typedef sparse_tag storage_category;
 841
 842         // Construction and destruction
 843         BOOST_UBLAS_INLINE
 844         compressed_vector ():
 845             vector_container<self_type> (),
 846             size_ (0), capacity_ (restrict_capacity (0)), filled_ (0),
 847             index_data_ (capacity_), value_data_ (capacity_) {
 848             storage_invariants ();
 849         }
 850         explicit BOOST_UBLAS_INLINE
 851         compressed_vector (size_type size, size_type non_zeros = 0):
 852             vector_container<self_type> (),
 853             size_ (size), capacity_ (restrict_capacity (non_zeros)), filled_ (0),
 854             index_data_ (capacity_), value_data_ (capacity_) {
 855         storage_invariants ();
 856         }
 857         BOOST_UBLAS_INLINE
 858         compressed_vector (const compressed_vector &v):
 859             vector_container<self_type> (),
 860             size_ (v.size_), capacity_ (v.capacity_), filled_ (v.filled_),
 861             index_data_ (v.index_data_), value_data_ (v.value_data_) {
 862             storage_invariants ();
 863         }
 864         template<class AE>
 865         BOOST_UBLAS_INLINE
 866         compressed_vector (const vector_expression<AE> &ae, size_type non_zeros = 0):
 867             vector_container<self_type> (),
 868             size_ (ae ().size ()), capacity_ (restrict_capacity (non_zeros)), filled_ (0),
 869             index_data_ (capacity_), value_data_ (capacity_) {
 870             storage_invariants ();
 871             vector_assign<scalar_assign> (*this, ae);
 872         }
 873
 874         // Accessors
 875         BOOST_UBLAS_INLINE
 876         size_type size () const {
 877             return size_;
 878         }
 879         BOOST_UBLAS_INLINE
 880         size_type nnz_capacity () const {
 881             return capacity_;
 882         }
 883         BOOST_UBLAS_INLINE
 884         size_type nnz () const {
 885             return filled_;
 886         }
 887
 888         // Storage accessors
 889         BOOST_UBLAS_INLINE
 890         static size_type index_base () {
 891             return IB;
 892         }
 893         BOOST_UBLAS_INLINE
 894         typename index_array_type::size_type filled () const {
 895             return filled_;
 896         }
 897         BOOST_UBLAS_INLINE
 898         const index_array_type &index_data () const {
 899             return index_data_;
 900         }
 901         BOOST_UBLAS_INLINE
 902         const value_array_type &value_data () const {
 903             return value_data_;
 904         }
 905         BOOST_UBLAS_INLINE
 906         void set_filled (const typename index_array_type::size_type & filled) {
 907             filled_ = filled;
 908             storage_invariants ();
 909         }
 910         BOOST_UBLAS_INLINE
 911         index_array_type &index_data () {
 912             return index_data_;
 913         }
 914         BOOST_UBLAS_INLINE
 915         value_array_type &value_data () {
 916             return value_data_;
 917         }
 918
 919         // Resizing
 920     private:
 921         BOOST_UBLAS_INLINE
 922         size_type restrict_capacity (size_type non_zeros) const {
 923             non_zeros = (std::max) (non_zeros, size_type (1));
 924             non_zeros = (std::min) (non_zeros, size_);
 925             return non_zeros;
 926         }
 927     public:
 928         BOOST_UBLAS_INLINE
 929         void resize (size_type size, bool preserve = true) {
 930             size_ = size;
 931             capacity_ = restrict_capacity (capacity_);
 932             if (preserve) {
 933                 index_data_. resize (capacity_, size_type ());
 934                 value_data_. resize (capacity_, value_type ());
 935                 filled_ = (std::min) (capacity_, filled_);
 936                 while ((filled_ > 0) && (zero_based(index_data_[filled_ - 1]) >= size)) {
 937                     --filled_;
 938                 }
 939             }
 940             else {
 941                 index_data_. resize (capacity_);
 942                 value_data_. resize (capacity_);
 943                 filled_ = 0;
 944             }
 945             storage_invariants ();
 946         }
 947
 948         // Reserving
 949         BOOST_UBLAS_INLINE
 950         void reserve (size_type non_zeros, bool preserve = true) {
 951             capacity_ = restrict_capacity (non_zeros);
 952             if (preserve) {
 953                 index_data_. resize (capacity_, size_type ());
 954                 value_data_. resize (capacity_, value_type ());
 955                 filled_ = (std::min) (capacity_, filled_);
 956             }
 957             else {
 958                 index_data_. resize (capacity_);
 959                 value_data_. resize (capacity_);
 960                 filled_ = 0;
 961             }
 962             storage_invariants ();
 963         }
 964
 965         // Element support
 966         BOOST_UBLAS_INLINE
 967         pointer find_element (size_type i) {
 968             return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i));
 969         }
 970         BOOST_UBLAS_INLINE
 971         const_pointer find_element (size_type i) const {
 972             const_subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
 973             if (it == index_data_.begin () + filled_ || *it != k_based (i))
 974                 return 0;
 975             return &value_data_ [it - index_data_.begin ()];
 976         }
 977
 978         // Element access
 979         BOOST_UBLAS_INLINE
 980         const_reference operator () (size_type i) const {
 981             BOOST_UBLAS_CHECK (i < size_, bad_index ());
 982             const_subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
 983             if (it == index_data_.begin () + filled_ || *it != k_based (i))
 984                 return zero_;
 985             return value_data_ [it - index_data_.begin ()];
 986         }
 987         BOOST_UBLAS_INLINE
 988         true_reference ref (size_type i) {
 989             BOOST_UBLAS_CHECK (i < size_, bad_index ());
 990             subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
 991             if (it == index_data_.begin () + filled_ || *it != k_based (i))
 992                 return insert_element (i, value_type/*zero*/());
 993             else
 994                 return value_data_ [it - index_data_.begin ()];
 995         }
 996         BOOST_UBLAS_INLINE
 997         reference operator () (size_type i) {
 998 #ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
 999             return ref (i) ;
1000 #else
1001             BOOST_UBLAS_CHECK (i < size_, bad_index ());
1002             return reference (*this, i);
1003 #endif
1004         }
1005
1006         BOOST_UBLAS_INLINE
1007         const_reference operator [] (size_type i) const {
1008             return (*this) (i);
1009         }
1010         BOOST_UBLAS_INLINE
1011         reference operator [] (size_type i) {
1012             return (*this) (i);
1013         }
1014
1015         // Element assignment
1016         BOOST_UBLAS_INLINE
1017         true_reference insert_element (size_type i, const_reference t) {
1018             BOOST_UBLAS_CHECK (!find_element (i), bad_index ());        // duplicate element
1019             if (filled_ >= capacity_)
1020                 reserve (2 * capacity_, true);
1021             subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1022             // ISSUE max_capacity limit due to difference_type
1023             typename std::iterator_traits<subiterator_type>::difference_type n = it - index_data_.begin ();
1024             BOOST_UBLAS_CHECK (filled_ == 0 || filled_ == typename index_array_type::size_type (n) || *it != k_based (i), internal_logic ());   // duplicate found by lower_bound
1025             ++ filled_;
1026             it = index_data_.begin () + n;
1027             std::copy_backward (it, index_data_.begin () + filled_ - 1, index_data_.begin () + filled_);
1028             *it = k_based (i);
1029             typename value_array_type::iterator itt (value_data_.begin () + n);
1030             std::copy_backward (itt, value_data_.begin () + filled_ - 1, value_data_.begin () + filled_);
1031             *itt = t;
1032             storage_invariants ();
1033             return *itt;
1034         }
1035         BOOST_UBLAS_INLINE
1036         void erase_element (size_type i) {
1037             subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1038             typename std::iterator_traits<subiterator_type>::difference_type  n = it - index_data_.begin ();
1039             if (filled_ > typename index_array_type::size_type (n) && *it == k_based (i)) {
1040                 std::copy (it + 1, index_data_.begin () + filled_, it);
1041                 typename value_array_type::iterator itt (value_data_.begin () + n);
1042                 std::copy (itt + 1, value_data_.begin () + filled_, itt);
1043                 -- filled_;
1044             }
1045             storage_invariants ();
1046         }
1047
1048         // Zeroing
1049         BOOST_UBLAS_INLINE
1050         void clear () {
1051             filled_ = 0;
1052             storage_invariants ();
1053         }
1054
1055         // Assignment
1056         BOOST_UBLAS_INLINE
1057         compressed_vector &operator = (const compressed_vector &v) {
1058             if (this != &v) {
1059                 size_ = v.size_;
1060                 capacity_ = v.capacity_;
1061                 filled_ = v.filled_;
1062                 index_data_ = v.index_data_;
1063                 value_data_ = v.value_data_;
1064             }
1065             storage_invariants ();
1066             return *this;
1067         }
1068         template<class C>          // Container assignment without temporary
1069         BOOST_UBLAS_INLINE
1070         compressed_vector &operator = (const vector_container<C> &v) {
1071             resize (v ().size (), false);
1072             assign (v);
1073             return *this;
1074         }
1075         BOOST_UBLAS_INLINE
1076         compressed_vector &assign_temporary (compressed_vector &v) {
1077             swap (v);
1078             return *this;
1079         }
1080         template<class AE>
1081         BOOST_UBLAS_INLINE
1082         compressed_vector &operator = (const vector_expression<AE> &ae) {
1083             self_type temporary (ae, capacity_);
1084             return assign_temporary (temporary);
1085         }
1086         template<class AE>
1087         BOOST_UBLAS_INLINE
1088         compressed_vector &assign (const vector_expression<AE> &ae) {
1089             vector_assign<scalar_assign> (*this, ae);
1090             return *this;
1091         }
1092
1093         // Computed assignment
1094         template<class AE>
1095         BOOST_UBLAS_INLINE
1096         compressed_vector &operator += (const vector_expression<AE> &ae) {
1097             self_type temporary (*this + ae, capacity_);
1098             return assign_temporary (temporary);
1099         }
1100         template<class C>          // Container assignment without temporary
1101         BOOST_UBLAS_INLINE
1102         compressed_vector &operator += (const vector_container<C> &v) {
1103             plus_assign (v);
1104             return *this;
1105         }
1106         template<class AE>
1107         BOOST_UBLAS_INLINE
1108         compressed_vector &plus_assign (const vector_expression<AE> &ae) {
1109             vector_assign<scalar_plus_assign> (*this, ae);
1110             return *this;
1111         }
1112         template<class AE>
1113         BOOST_UBLAS_INLINE
1114         compressed_vector &operator -= (const vector_expression<AE> &ae) {
1115             self_type temporary (*this - ae, capacity_);
1116             return assign_temporary (temporary);
1117         }
1118         template<class C>          // Container assignment without temporary
1119         BOOST_UBLAS_INLINE
1120         compressed_vector &operator -= (const vector_container<C> &v) {
1121             minus_assign (v);
1122             return *this;
1123         }
1124         template<class AE>
1125         BOOST_UBLAS_INLINE
1126         compressed_vector &minus_assign (const vector_expression<AE> &ae) {
1127             vector_assign<scalar_minus_assign> (*this, ae);
1128             return *this;
1129         }
1130         template<class AT>
1131         BOOST_UBLAS_INLINE
1132         compressed_vector &operator *= (const AT &at) {
1133             vector_assign_scalar<scalar_multiplies_assign> (*this, at);
1134             return *this;
1135         }
1136         template<class AT>
1137         BOOST_UBLAS_INLINE
1138         compressed_vector &operator /= (const AT &at) {
1139             vector_assign_scalar<scalar_divides_assign> (*this, at);
1140             return *this;
1141         }
1142
1143         // Swapping
1144         BOOST_UBLAS_INLINE
1145         void swap (compressed_vector &v) {
1146             if (this != &v) {
1147                 std::swap (size_, v.size_);
1148                 std::swap (capacity_, v.capacity_);
1149                 std::swap (filled_, v.filled_);
1150                 index_data_.swap (v.index_data_);
1151                 value_data_.swap (v.value_data_);
1152             }
1153             storage_invariants ();
1154         }
1155         BOOST_UBLAS_INLINE
1156         friend void swap (compressed_vector &v1, compressed_vector &v2) {
1157             v1.swap (v2);
1158         }
1159
1160         // Back element insertion and erasure
1161         BOOST_UBLAS_INLINE
1162         void push_back (size_type i, const_reference t) {
1163             BOOST_UBLAS_CHECK (filled_ == 0 || index_data_ [filled_ - 1] < k_based (i), external_logic ());
1164             if (filled_ >= capacity_)
1165                 reserve (2 * capacity_, true);
1166             BOOST_UBLAS_CHECK (filled_ < capacity_, internal_logic ());
1167             index_data_ [filled_] = k_based (i);
1168             value_data_ [filled_] = t;
1169             ++ filled_;
1170             storage_invariants ();
1171         }
1172         BOOST_UBLAS_INLINE
1173         void pop_back () {
1174             BOOST_UBLAS_CHECK (filled_ > 0, external_logic ());
1175             -- filled_;
1176             storage_invariants ();
1177         }
1178
1179         // Iterator types
1180     private:
1181         // Use index array iterator
1182         typedef typename IA::const_iterator const_subiterator_type;
1183         typedef typename IA::iterator subiterator_type;
1184
1185         BOOST_UBLAS_INLINE
1186         true_reference at_element (size_type i) {
1187             BOOST_UBLAS_CHECK (i < size_, bad_index ());
1188             subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1189             BOOST_UBLAS_CHECK (it != index_data_.begin () + filled_ && *it == k_based (i), bad_index ());
1190             return value_data_ [it - index_data_.begin ()];
1191         }
1192
1193     public:
1194         class const_iterator;
1195         class iterator;
1196
1197         // Element lookup
1198         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1199         const_iterator find (size_type i) const {
1200             return const_iterator (*this, detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1201         }
1202         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1203         iterator find (size_type i) {
1204             return iterator (*this, detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1205         }
1206
1207
1208         class const_iterator:
1209             public container_const_reference<compressed_vector>,
1210             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
1211                                                const_iterator, value_type> {
1212         public:
1213             typedef typename compressed_vector::value_type value_type;
1214             typedef typename compressed_vector::difference_type difference_type;
1215             typedef typename compressed_vector::const_reference reference;
1216             typedef const typename compressed_vector::pointer pointer;
1217
1218             // Construction and destruction
1219             BOOST_UBLAS_INLINE
1220             const_iterator ():
1221                 container_const_reference<self_type> (), it_ () {}
1222             BOOST_UBLAS_INLINE
1223             const_iterator (const self_type &v, const const_subiterator_type &it):
1224                 container_const_reference<self_type> (v), it_ (it) {}
1225             BOOST_UBLAS_INLINE
1226             const_iterator (const typename self_type::iterator &it):  // ISSUE self_type:: stops VC8 using std::iterator here
1227                 container_const_reference<self_type> (it ()), it_ (it.it_) {}
1228
1229             // Arithmetic
1230             BOOST_UBLAS_INLINE
1231             const_iterator &operator ++ () {
1232                 ++ it_;
1233                 return *this;
1234             }
1235             BOOST_UBLAS_INLINE
1236             const_iterator &operator -- () {
1237                 -- it_;
1238                 return *this;
1239             }
1240
1241             // Dereference
1242             BOOST_UBLAS_INLINE
1243             const_reference operator * () const {
1244                 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
1245                 return (*this) ().value_data_ [it_ - (*this) ().index_data_.begin ()];
1246             }
1247
1248             // Index
1249             BOOST_UBLAS_INLINE
1250             size_type index () const {
1251                 BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
1252                 BOOST_UBLAS_CHECK ((*this) ().zero_based (*it_) < (*this) ().size (), bad_index ());
1253                 return (*this) ().zero_based (*it_);
1254             }
1255
1256             // Assignment
1257             BOOST_UBLAS_INLINE
1258             const_iterator &operator = (const const_iterator &it) {
1259                 container_const_reference<self_type>::assign (&it ());
1260                 it_ = it.it_;
1261                 return *this;
1262             }
1263
1264             // Comparison
1265             BOOST_UBLAS_INLINE
1266             bool operator == (const const_iterator &it) const {
1267                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1268                 return it_ == it.it_;
1269             }
1270
1271         private:
1272             const_subiterator_type it_;
1273         };
1274
1275         BOOST_UBLAS_INLINE
1276         const_iterator begin () const {
1277             return find (0);
1278         }
1279         BOOST_UBLAS_INLINE
1280         const_iterator cbegin () const {
1281             return begin ();
1282         }
1283         BOOST_UBLAS_INLINE
1284         const_iterator end () const {
1285             return find (size_);
1286         }
1287         BOOST_UBLAS_INLINE
1288         const_iterator cend () const {
1289             return end ();
1290         }
1291
1292         class iterator:
1293             public container_reference<compressed_vector>,
1294             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
1295                                                iterator, value_type> {
1296         public:
1297             typedef typename compressed_vector::value_type value_type;
1298             typedef typename compressed_vector::difference_type difference_type;
1299             typedef typename compressed_vector::true_reference reference;
1300             typedef typename compressed_vector::pointer pointer;
1301
1302             // Construction and destruction
1303             BOOST_UBLAS_INLINE
1304             iterator ():
1305                 container_reference<self_type> (), it_ () {}
1306             BOOST_UBLAS_INLINE
1307             iterator (self_type &v, const subiterator_type &it):
1308                 container_reference<self_type> (v), it_ (it) {}
1309
1310             // Arithmetic
1311             BOOST_UBLAS_INLINE
1312             iterator &operator ++ () {
1313                 ++ it_;
1314                 return *this;
1315             }
1316             BOOST_UBLAS_INLINE
1317             iterator &operator -- () {
1318                 -- it_;
1319                 return *this;
1320             }
1321
1322             // Dereference
1323             BOOST_UBLAS_INLINE
1324             reference operator * () const {
1325                 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
1326                 return (*this) ().value_data_ [it_ - (*this) ().index_data_.begin ()];
1327             }
1328
1329             // Index
1330             BOOST_UBLAS_INLINE
1331             size_type index () const {
1332                 BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
1333                 BOOST_UBLAS_CHECK ((*this) ().zero_based (*it_) < (*this) ().size (), bad_index ());
1334                 return (*this) ().zero_based (*it_);
1335             }
1336
1337             // Assignment
1338             BOOST_UBLAS_INLINE
1339             iterator &operator = (const iterator &it) {
1340                 container_reference<self_type>::assign (&it ());
1341                 it_ = it.it_;
1342                 return *this;
1343             }
1344
1345             // Comparison
1346             BOOST_UBLAS_INLINE
1347             bool operator == (const iterator &it) const {
1348                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1349                 return it_ == it.it_;
1350             }
1351
1352         private:
1353             subiterator_type it_;
1354
1355             friend class const_iterator;
1356         };
1357
1358         BOOST_UBLAS_INLINE
1359         iterator begin () {
1360             return find (0);
1361         }
1362         BOOST_UBLAS_INLINE
1363         iterator end () {
1364             return find (size_);
1365         }
1366
1367         // Reverse iterator
1368         typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
1369         typedef reverse_iterator_base<iterator> reverse_iterator;
1370
1371         BOOST_UBLAS_INLINE
1372         const_reverse_iterator rbegin () const {
1373             return const_reverse_iterator (end ());
1374         }
1375         BOOST_UBLAS_INLINE
1376         const_reverse_iterator crbegin () const {
1377             return rbegin ();
1378         }
1379         BOOST_UBLAS_INLINE
1380         const_reverse_iterator rend () const {
1381             return const_reverse_iterator (begin ());
1382         }
1383         BOOST_UBLAS_INLINE
1384         const_reverse_iterator crend () const {
1385             return rend ();
1386         }
1387         BOOST_UBLAS_INLINE
1388         reverse_iterator rbegin () {
1389             return reverse_iterator (end ());
1390         }
1391         BOOST_UBLAS_INLINE
1392         reverse_iterator rend () {
1393             return reverse_iterator (begin ());
1394         }
1395
1396          // Serialization
1397         template<class Archive>
1398         void serialize(Archive & ar, const unsigned int /* file_version */){
1399             serialization::collection_size_type s (size_);
1400             ar & serialization::make_nvp("size",s);
1401             if (Archive::is_loading::value) {
1402                 size_ = s;
1403             }
1404             // ISSUE: filled may be much less than capacity
1405             // ISSUE: index_data_ and value_data_ are undefined between filled and capacity (trouble with 'nan'-values)
1406             ar & serialization::make_nvp("capacity", capacity_);
1407             ar & serialization::make_nvp("filled", filled_);
1408             ar & serialization::make_nvp("index_data", index_data_);
1409             ar & serialization::make_nvp("value_data", value_data_);
1410             storage_invariants();
1411         }
1412
1413     private:
1414         void storage_invariants () const
1415         {
1416             BOOST_UBLAS_CHECK (capacity_ == index_data_.size (), internal_logic ());
1417             BOOST_UBLAS_CHECK (capacity_ == value_data_.size (), internal_logic ());
1418             BOOST_UBLAS_CHECK (filled_ <= capacity_, internal_logic ());
1419             BOOST_UBLAS_CHECK ((0 == filled_) || (zero_based(index_data_[filled_ - 1]) < size_), internal_logic ());
1420         }
1421
1422         size_type size_;
1423         typename index_array_type::size_type capacity_;
1424         typename index_array_type::size_type filled_;
1425         index_array_type index_data_;
1426         value_array_type value_data_;
1427         static const value_type zero_;
1428
1429         BOOST_UBLAS_INLINE
1430         static size_type zero_based (size_type k_based_index) {
1431             return k_based_index - IB;
1432         }
1433         BOOST_UBLAS_INLINE
1434         static size_type k_based (size_type zero_based_index) {
1435             return zero_based_index + IB;
1436         }
1437
1438         friend class iterator;
1439         friend class const_iterator;
1440     };
1441
1442     template<class T, std::size_t IB, class IA, class TA>
1443     const typename compressed_vector<T, IB, IA, TA>::value_type compressed_vector<T, IB, IA, TA>::zero_ = value_type/*zero*/();
1444
1445     // Thanks to Kresimir Fresl for extending this to cover different index bases.
1446
1447     /** \brief Coordimate array based sparse vector
1448      *
1449      * a sparse vector of values of type \c T of variable size. The non zero values are stored
1450      * as two seperate arrays: an index array and a value array. The arrays may be out of order
1451      * with multiple entries for each vector element. If there are multiple values for the same
1452      * index the sum of these values is the real value. It is way more efficient for inserting values
1453      * than a \c compressed_vector but less memory efficient. Also linearly parsing a vector can
1454      * be longer in specific cases than a \c compressed_vector.
1455      *
1456      * For a n-dimensional sorted coordinate vector and \f$ 0 \leq i < n\f$ the non-zero elements
1457      * \f$v_i\f$ are mapped to consecutive elements of the index and value container, i.e. for
1458      * 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$.
1459      *
1460      * Supported parameters for the adapted array (indices and values) are \c unbounded_array<> ,
1461      * \c bounded_array<> and \c std::vector<>.
1462      *
1463      * \tparam T the type of object stored in the vector (like double, float, complex, etc...)
1464      * \tparam IB the index base of the compressed vector. Default is 0. Other supported value is 1
1465      * \tparam IA the type of adapted array for indices. Default is \c unbounded_array<std::size_t>
1466      * \tparam TA the type of adapted array for values. Default is unbounded_array<T>
1467      */
1468     template<class T, std::size_t IB, class IA, class TA>
1469     class coordinate_vector:
1470         public vector_container<coordinate_vector<T, IB, IA, TA> > {
1471
1472         typedef T &true_reference;
1473         typedef T *pointer;
1474         typedef const T *const_pointer;
1475         typedef coordinate_vector<T, IB, IA, TA> self_type;
1476     public:
1477 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
1478         using vector_container<self_type>::operator ();
1479 #endif
1480         // ISSUE require type consistency check
1481         // is_convertable (IA::size_type, TA::size_type)
1482         typedef typename IA::value_type size_type;
1483         typedef typename IA::difference_type difference_type;
1484         typedef T value_type;
1485         typedef const T &const_reference;
1486 #ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
1487         typedef T &reference;
1488 #else
1489         typedef sparse_vector_element<self_type> reference;
1490 #endif
1491         typedef IA index_array_type;
1492         typedef TA value_array_type;
1493         typedef const vector_reference<const self_type> const_closure_type;
1494         typedef vector_reference<self_type> closure_type;
1495         typedef self_type vector_temporary_type;
1496         typedef sparse_tag storage_category;
1497
1498         // Construction and destruction
1499         BOOST_UBLAS_INLINE
1500         coordinate_vector ():
1501             vector_container<self_type> (),
1502             size_ (0), capacity_ (restrict_capacity (0)),
1503             filled_ (0), sorted_filled_ (filled_), sorted_ (true),
1504             index_data_ (capacity_), value_data_ (capacity_) {
1505             storage_invariants ();
1506         }
1507         explicit BOOST_UBLAS_INLINE
1508         coordinate_vector (size_type size, size_type non_zeros = 0):
1509             vector_container<self_type> (),
1510             size_ (size), capacity_ (restrict_capacity (non_zeros)),
1511             filled_ (0), sorted_filled_ (filled_), sorted_ (true),
1512             index_data_ (capacity_), value_data_ (capacity_) {
1513             storage_invariants ();
1514         }
1515         BOOST_UBLAS_INLINE
1516         coordinate_vector (const coordinate_vector &v):
1517             vector_container<self_type> (),
1518             size_ (v.size_), capacity_ (v.capacity_),
1519             filled_ (v.filled_), sorted_filled_ (v.sorted_filled_), sorted_ (v.sorted_),
1520             index_data_ (v.index_data_), value_data_ (v.value_data_) {
1521             storage_invariants ();
1522         }
1523         template<class AE>
1524         BOOST_UBLAS_INLINE
1525         coordinate_vector (const vector_expression<AE> &ae, size_type non_zeros = 0):
1526             vector_container<self_type> (),
1527             size_ (ae ().size ()), capacity_ (restrict_capacity (non_zeros)),
1528             filled_ (0), sorted_filled_ (filled_), sorted_ (true),
1529             index_data_ (capacity_), value_data_ (capacity_) {
1530             storage_invariants ();
1531             vector_assign<scalar_assign> (*this, ae);
1532         }
1533
1534         // Accessors
1535         BOOST_UBLAS_INLINE
1536         size_type size () const {
1537             return size_;
1538         }
1539         BOOST_UBLAS_INLINE
1540         size_type nnz_capacity () const {
1541             return capacity_;
1542         }
1543         BOOST_UBLAS_INLINE
1544         size_type nnz () const {
1545             return filled_;
1546         }
1547
1548         // Storage accessors
1549         BOOST_UBLAS_INLINE
1550         static size_type index_base () {
1551             return IB;
1552         }
1553         BOOST_UBLAS_INLINE
1554         typename index_array_type::size_type filled () const {
1555             return filled_;
1556         }
1557         BOOST_UBLAS_INLINE
1558         const index_array_type &index_data () const {
1559             return index_data_;
1560         }
1561         BOOST_UBLAS_INLINE
1562         const value_array_type &value_data () const {
1563             return value_data_;
1564         }
1565         BOOST_UBLAS_INLINE
1566         void set_filled (const typename index_array_type::size_type &sorted, const typename index_array_type::size_type &filled) {
1567             sorted_filled_ = sorted;
1568             filled_ = filled;
1569             storage_invariants ();
1570         }
1571         BOOST_UBLAS_INLINE
1572         index_array_type &index_data () {
1573             return index_data_;
1574         }
1575         BOOST_UBLAS_INLINE
1576         value_array_type &value_data () {
1577             return value_data_;
1578         }
1579
1580         // Resizing
1581     private:
1582         BOOST_UBLAS_INLINE
1583         size_type restrict_capacity (size_type non_zeros) const {
1584              // minimum non_zeros
1585              non_zeros = (std::max) (non_zeros, size_type (1));
1586              // ISSUE no maximum as coordinate may contain inserted duplicates
1587              return non_zeros;
1588         }
1589     public:
1590         BOOST_UBLAS_INLINE
1591         void resize (size_type size, bool preserve = true) {
1592             if (preserve)
1593                 sort ();    // remove duplicate elements.
1594             size_ = size;
1595             capacity_ = restrict_capacity (capacity_);
1596             if (preserve) {
1597                 index_data_. resize (capacity_, size_type ());
1598                 value_data_. resize (capacity_, value_type ());
1599                 filled_ = (std::min) (capacity_, filled_);
1600                 while ((filled_ > 0) && (zero_based(index_data_[filled_ - 1]) >= size)) {
1601                     --filled_;
1602                 }
1603             }
1604             else {
1605                 index_data_. resize (capacity_);
1606                 value_data_. resize (capacity_);
1607                 filled_ = 0;
1608             }
1609             sorted_filled_ = filled_;
1610             storage_invariants ();
1611         }
1612         // Reserving
1613         BOOST_UBLAS_INLINE
1614         void reserve (size_type non_zeros, bool preserve = true) {
1615             if (preserve)
1616                 sort ();    // remove duplicate elements.
1617             capacity_ = restrict_capacity (non_zeros);
1618             if (preserve) {
1619                 index_data_. resize (capacity_, size_type ());
1620                 value_data_. resize (capacity_, value_type ());
1621                 filled_ = (std::min) (capacity_, filled_);
1622                 }
1623             else {
1624                 index_data_. resize (capacity_);
1625                 value_data_. resize (capacity_);
1626                 filled_ = 0;
1627             }
1628             sorted_filled_ = filled_;
1629             storage_invariants ();
1630         }
1631
1632         // Element support
1633         BOOST_UBLAS_INLINE
1634         pointer find_element (size_type i) {
1635             return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i));
1636         }
1637         BOOST_UBLAS_INLINE
1638         const_pointer find_element (size_type i) const {
1639             sort ();
1640             const_subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1641             if (it == index_data_.begin () + filled_ || *it != k_based (i))
1642                 return 0;
1643             return &value_data_ [it - index_data_.begin ()];
1644         }
1645
1646         // Element access
1647         BOOST_UBLAS_INLINE
1648         const_reference operator () (size_type i) const {
1649             BOOST_UBLAS_CHECK (i < size_, bad_index ());
1650             sort ();
1651             const_subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1652             if (it == index_data_.begin () + filled_ || *it != k_based (i))
1653                 return zero_;
1654             return value_data_ [it - index_data_.begin ()];
1655         }
1656         BOOST_UBLAS_INLINE
1657         true_reference ref (size_type i) {
1658             BOOST_UBLAS_CHECK (i < size_, bad_index ());
1659             sort ();
1660             subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1661             if (it == index_data_.begin () + filled_ || *it != k_based (i))
1662                 return insert_element (i, value_type/*zero*/());
1663             else
1664                 return value_data_ [it - index_data_.begin ()];
1665         }
1666         BOOST_UBLAS_INLINE
1667         reference operator () (size_type i) {
1668 #ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
1669             return ref (i);
1670 #else
1671             BOOST_UBLAS_CHECK (i < size_, bad_index ());
1672             return reference (*this, i);
1673 #endif
1674         }
1675
1676         BOOST_UBLAS_INLINE
1677         const_reference operator [] (size_type i) const {
1678             return (*this) (i);
1679         }
1680         BOOST_UBLAS_INLINE
1681         reference operator [] (size_type i) {
1682             return (*this) (i);
1683         }
1684
1685         // Element assignment
1686         BOOST_UBLAS_INLINE
1687         void append_element (size_type i, const_reference t) {
1688             if (filled_ >= capacity_)
1689                 reserve (2 * filled_, true);
1690             BOOST_UBLAS_CHECK (filled_ < capacity_, internal_logic ());
1691             index_data_ [filled_] = k_based (i);
1692             value_data_ [filled_] = t;
1693             ++ filled_;
1694             sorted_ = false;
1695             storage_invariants ();
1696         }
1697         BOOST_UBLAS_INLINE
1698         true_reference insert_element (size_type i, const_reference t) {
1699             BOOST_UBLAS_CHECK (!find_element (i), bad_index ());        // duplicate element
1700             append_element (i, t);
1701             return value_data_ [filled_ - 1];
1702         }
1703         BOOST_UBLAS_INLINE
1704         void erase_element (size_type i) {
1705             sort ();
1706             subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1707             typename std::iterator_traits<subiterator_type>::difference_type n = it - index_data_.begin ();
1708             if (filled_ > typename index_array_type::size_type (n) && *it == k_based (i)) {
1709                 std::copy (it + 1, index_data_.begin () + filled_, it);
1710                 typename value_array_type::iterator itt (value_data_.begin () + n);
1711                 std::copy (itt + 1, value_data_.begin () + filled_, itt);
1712                 -- filled_;
1713                 sorted_filled_ = filled_;
1714             }
1715             storage_invariants ();
1716         }
1717
1718         // Zeroing
1719         BOOST_UBLAS_INLINE
1720         void clear () {
1721             filled_ = 0;
1722             sorted_filled_ = filled_;
1723             sorted_ = true;
1724             storage_invariants ();
1725         }
1726
1727         // Assignment
1728         BOOST_UBLAS_INLINE
1729         coordinate_vector &operator = (const coordinate_vector &v) {
1730             if (this != &v) {
1731                 size_ = v.size_;
1732                 capacity_ = v.capacity_;
1733                 filled_ = v.filled_;
1734                 sorted_filled_ = v.sorted_filled_;
1735                 sorted_ = v.sorted_;
1736                 index_data_ = v.index_data_;
1737                 value_data_ = v.value_data_;
1738             }
1739             storage_invariants ();
1740             return *this;
1741         }
1742         template<class C>          // Container assignment without temporary
1743         BOOST_UBLAS_INLINE
1744         coordinate_vector &operator = (const vector_container<C> &v) {
1745             resize (v ().size (), false);
1746             assign (v);
1747             return *this;
1748         }
1749         BOOST_UBLAS_INLINE
1750         coordinate_vector &assign_temporary (coordinate_vector &v) {
1751             swap (v);
1752             return *this;
1753         }
1754         template<class AE>
1755         BOOST_UBLAS_INLINE
1756         coordinate_vector &operator = (const vector_expression<AE> &ae) {
1757             self_type temporary (ae, capacity_);
1758             return assign_temporary (temporary);
1759         }
1760         template<class AE>
1761         BOOST_UBLAS_INLINE
1762         coordinate_vector &assign (const vector_expression<AE> &ae) {
1763             vector_assign<scalar_assign> (*this, ae);
1764             return *this;
1765         }
1766
1767         // Computed assignment
1768         template<class AE>
1769         BOOST_UBLAS_INLINE
1770         coordinate_vector &operator += (const vector_expression<AE> &ae) {
1771             self_type temporary (*this + ae, capacity_);
1772             return assign_temporary (temporary);
1773         }
1774         template<class C>          // Container assignment without temporary
1775         BOOST_UBLAS_INLINE
1776         coordinate_vector &operator += (const vector_container<C> &v) {
1777             plus_assign (v);
1778             return *this;
1779         }
1780         template<class AE>
1781         BOOST_UBLAS_INLINE
1782         coordinate_vector &plus_assign (const vector_expression<AE> &ae) {
1783             vector_assign<scalar_plus_assign> (*this, ae);
1784             return *this;
1785         }
1786         template<class AE>
1787         BOOST_UBLAS_INLINE
1788         coordinate_vector &operator -= (const vector_expression<AE> &ae) {
1789             self_type temporary (*this - ae, capacity_);
1790             return assign_temporary (temporary);
1791         }
1792         template<class C>          // Container assignment without temporary
1793         BOOST_UBLAS_INLINE
1794         coordinate_vector &operator -= (const vector_container<C> &v) {
1795             minus_assign (v);
1796             return *this;
1797         }
1798         template<class AE>
1799         BOOST_UBLAS_INLINE
1800         coordinate_vector &minus_assign (const vector_expression<AE> &ae) {
1801             vector_assign<scalar_minus_assign> (*this, ae);
1802             return *this;
1803         }
1804         template<class AT>
1805         BOOST_UBLAS_INLINE
1806         coordinate_vector &operator *= (const AT &at) {
1807             vector_assign_scalar<scalar_multiplies_assign> (*this, at);
1808             return *this;
1809         }
1810         template<class AT>
1811         BOOST_UBLAS_INLINE
1812         coordinate_vector &operator /= (const AT &at) {
1813             vector_assign_scalar<scalar_divides_assign> (*this, at);
1814             return *this;
1815         }
1816
1817         // Swapping
1818         BOOST_UBLAS_INLINE
1819         void swap (coordinate_vector &v) {
1820             if (this != &v) {
1821                 std::swap (size_, v.size_);
1822                 std::swap (capacity_, v.capacity_);
1823                 std::swap (filled_, v.filled_);
1824                 std::swap (sorted_filled_, v.sorted_filled_);
1825                 std::swap (sorted_, v.sorted_);
1826                 index_data_.swap (v.index_data_);
1827                 value_data_.swap (v.value_data_);
1828             }
1829             storage_invariants ();
1830         }
1831         BOOST_UBLAS_INLINE
1832         friend void swap (coordinate_vector &v1, coordinate_vector &v2) {
1833             v1.swap (v2);
1834         }
1835
1836         // replacement if STL lower bound algorithm for use of inplace_merge
1837         size_type lower_bound (size_type beg, size_type end, size_type target) const {
1838             while (end > beg) {
1839                 size_type mid = (beg + end) / 2;
1840                 if (index_data_[mid] < index_data_[target]) {
1841                     beg = mid + 1;
1842                 } else {
1843                     end = mid;
1844                 }
1845             }
1846             return beg;
1847         }
1848
1849         // specialized replacement of STL inplace_merge to avoid compilation
1850         // problems with respect to the array_triple iterator
1851         void inplace_merge (size_type beg, size_type mid, size_type end) const {
1852             size_type len_lef = mid - beg;
1853             size_type len_rig = end - mid;
1854
1855             if (len_lef == 1 && len_rig == 1) {
1856                 if (index_data_[mid] < index_data_[beg]) {
1857                     std::swap(index_data_[beg], index_data_[mid]);
1858                     std::swap(value_data_[beg], value_data_[mid]);
1859                 }
1860             } else if (len_lef > 0 && len_rig > 0) {
1861                 size_type lef_mid, rig_mid;
1862                 if (len_lef >= len_rig) {
1863                     lef_mid = (beg + mid) / 2;
1864                     rig_mid = lower_bound(mid, end, lef_mid);
1865                 } else {
1866                     rig_mid = (mid + end) / 2;
1867                     lef_mid = lower_bound(beg, mid, rig_mid);
1868                 }
1869                 std::rotate(&index_data_[0] + lef_mid, &index_data_[0] + mid, &index_data_[0] + rig_mid);
1870                 std::rotate(&value_data_[0] + lef_mid, &value_data_[0] + mid, &value_data_[0] + rig_mid);
1871
1872                 size_type new_mid = lef_mid + rig_mid - mid;
1873                 inplace_merge(beg, lef_mid, new_mid);
1874                 inplace_merge(new_mid, rig_mid, end);
1875             }
1876         }
1877
1878         // Sorting and summation of duplicates
1879         BOOST_UBLAS_INLINE
1880         void sort () const {
1881             if (! sorted_ && filled_ > 0) {
1882                 typedef index_pair_array<index_array_type, value_array_type> array_pair;
1883                 array_pair ipa (filled_, index_data_, value_data_);
1884 #ifndef BOOST_UBLAS_COO_ALWAYS_DO_FULL_SORT
1885                 const typename array_pair::iterator iunsorted = ipa.begin () + sorted_filled_;
1886                 // sort new elements and merge
1887                 std::sort (iunsorted, ipa.end ());
1888                 inplace_merge(0, sorted_filled_, filled_);
1889 #else
1890                 const typename array_pair::iterator iunsorted = ipa.begin ();
1891                 std::sort (iunsorted, ipa.end ());
1892 #endif
1893
1894                 // sum duplicates with += and remove
1895                 size_type filled = 0;
1896                 for (size_type i = 1; i < filled_; ++ i) {
1897                     if (index_data_ [filled] != index_data_ [i]) {
1898                         ++ filled;
1899                         if (filled != i) {
1900                             index_data_ [filled] = index_data_ [i];
1901                             value_data_ [filled] = value_data_ [i];
1902                         }
1903                     } else {
1904                         value_data_ [filled] += value_data_ [i];
1905                     }
1906                 }
1907                 filled_ = filled + 1;
1908                 sorted_filled_ = filled_;
1909                 sorted_ = true;
1910                 storage_invariants ();
1911             }
1912         }
1913
1914         // Back element insertion and erasure
1915         BOOST_UBLAS_INLINE
1916         void push_back (size_type i, const_reference t) {
1917             // must maintain sort order
1918             BOOST_UBLAS_CHECK (sorted_ && (filled_ == 0 || index_data_ [filled_ - 1] < k_based (i)), external_logic ());
1919             if (filled_ >= capacity_)
1920                 reserve (2 * filled_, true);
1921             BOOST_UBLAS_CHECK (filled_ < capacity_, internal_logic ());
1922             index_data_ [filled_] = k_based (i);
1923             value_data_ [filled_] = t;
1924             ++ filled_;
1925             sorted_filled_ = filled_;
1926             storage_invariants ();
1927         }
1928         BOOST_UBLAS_INLINE
1929         void pop_back () {
1930             // ISSUE invariants could be simpilfied if sorted required as precondition
1931             BOOST_UBLAS_CHECK (filled_ > 0, external_logic ());
1932             -- filled_;
1933             sorted_filled_ = (std::min) (sorted_filled_, filled_);
1934             sorted_ = sorted_filled_ = filled_;
1935             storage_invariants ();
1936         }
1937
1938         // Iterator types
1939     private:
1940         // Use index array iterator
1941         typedef typename IA::const_iterator const_subiterator_type;
1942         typedef typename IA::iterator subiterator_type;
1943
1944         BOOST_UBLAS_INLINE
1945         true_reference at_element (size_type i) {
1946             BOOST_UBLAS_CHECK (i < size_, bad_index ());
1947             sort ();
1948             subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1949             BOOST_UBLAS_CHECK (it != index_data_.begin () + filled_ && *it == k_based (i), bad_index ());
1950             return value_data_ [it - index_data_.begin ()];
1951         }
1952
1953     public:
1954         class const_iterator;
1955         class iterator;
1956
1957         // Element lookup
1958         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1959         const_iterator find (size_type i) const {
1960             sort ();
1961             return const_iterator (*this, detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1962         }
1963         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1964         iterator find (size_type i) {
1965             sort ();
1966             return iterator (*this, detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1967         }
1968
1969
1970         class const_iterator:
1971             public container_const_reference<coordinate_vector>,
1972             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
1973                                                const_iterator, value_type> {
1974         public:
1975             typedef typename coordinate_vector::value_type value_type;
1976             typedef typename coordinate_vector::difference_type difference_type;
1977             typedef typename coordinate_vector::const_reference reference;
1978             typedef const typename coordinate_vector::pointer pointer;
1979
1980             // Construction and destruction
1981             BOOST_UBLAS_INLINE
1982             const_iterator ():
1983                 container_const_reference<self_type> (), it_ () {}
1984             BOOST_UBLAS_INLINE
1985             const_iterator (const self_type &v, const const_subiterator_type &it):
1986                 container_const_reference<self_type> (v), it_ (it) {}
1987             BOOST_UBLAS_INLINE
1988             const_iterator (const typename self_type::iterator &it):  // ISSUE self_type:: stops VC8 using std::iterator here
1989                 container_const_reference<self_type> (it ()), it_ (it.it_) {}
1990
1991             // Arithmetic
1992             BOOST_UBLAS_INLINE
1993             const_iterator &operator ++ () {
1994                 ++ it_;
1995                 return *this;
1996             }
1997             BOOST_UBLAS_INLINE
1998             const_iterator &operator -- () {
1999                 -- it_;
2000                 return *this;
2001             }
2002
2003             // Dereference
2004             BOOST_UBLAS_INLINE
2005             const_reference operator * () const {
2006                 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
2007                 return (*this) ().value_data_ [it_ - (*this) ().index_data_.begin ()];
2008             }
2009
2010             // Index
2011             BOOST_UBLAS_INLINE
2012             size_type index () const {
2013                 BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
2014                 BOOST_UBLAS_CHECK ((*this) ().zero_based (*it_) < (*this) ().size (), bad_index ());
2015                 return (*this) ().zero_based (*it_);
2016             }
2017
2018             // Assignment
2019             BOOST_UBLAS_INLINE
2020             const_iterator &operator = (const const_iterator &it) {
2021                 container_const_reference<self_type>::assign (&it ());
2022                 it_ = it.it_;
2023                 return *this;
2024             }
2025
2026             // Comparison
2027             BOOST_UBLAS_INLINE
2028             bool operator == (const const_iterator &it) const {
2029                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2030                 return it_ == it.it_;
2031             }
2032
2033         private:
2034             const_subiterator_type it_;
2035         };
2036
2037         BOOST_UBLAS_INLINE
2038         const_iterator begin () const {
2039             return find (0);
2040         }
2041         BOOST_UBLAS_INLINE
2042         const_iterator cbegin () const {
2043             return begin();
2044         }
2045         BOOST_UBLAS_INLINE
2046         const_iterator end () const {
2047             return find (size_);
2048         }
2049         BOOST_UBLAS_INLINE
2050         const_iterator cend () const {
2051             return end();
2052         }
2053
2054         class iterator:
2055             public container_reference<coordinate_vector>,
2056             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
2057                                                iterator, value_type> {
2058         public:
2059             typedef typename coordinate_vector::value_type value_type;
2060             typedef typename coordinate_vector::difference_type difference_type;
2061             typedef typename coordinate_vector::true_reference reference;
2062             typedef typename coordinate_vector::pointer pointer;
2063
2064             // Construction and destruction
2065             BOOST_UBLAS_INLINE
2066             iterator ():
2067                 container_reference<self_type> (), it_ () {}
2068             BOOST_UBLAS_INLINE
2069             iterator (self_type &v, const subiterator_type &it):
2070                 container_reference<self_type> (v), it_ (it) {}
2071
2072             // Arithmetic
2073             BOOST_UBLAS_INLINE
2074             iterator &operator ++ () {
2075                 ++ it_;
2076                 return *this;
2077             }
2078             BOOST_UBLAS_INLINE
2079             iterator &operator -- () {
2080                 -- it_;
2081                 return *this;
2082             }
2083
2084             // Dereference
2085             BOOST_UBLAS_INLINE
2086             reference operator * () const {
2087                 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
2088                 return (*this) ().value_data_ [it_ - (*this) ().index_data_.begin ()];
2089             }
2090
2091             // Index
2092             BOOST_UBLAS_INLINE
2093             size_type index () const {
2094                 BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
2095                 BOOST_UBLAS_CHECK ((*this) ().zero_based (*it_) < (*this) ().size (), bad_index ());
2096                 return (*this) ().zero_based (*it_);
2097             }
2098
2099             // Assignment
2100             BOOST_UBLAS_INLINE
2101             iterator &operator = (const iterator &it) {
2102                 container_reference<self_type>::assign (&it ());
2103                 it_ = it.it_;
2104                 return *this;
2105             }
2106
2107             // Comparison
2108             BOOST_UBLAS_INLINE
2109             bool operator == (const iterator &it) const {
2110                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2111                 return it_ == it.it_;
2112             }
2113
2114         private:
2115             subiterator_type it_;
2116
2117             friend class const_iterator;
2118         };
2119
2120         BOOST_UBLAS_INLINE
2121         iterator begin () {
2122             return find (0);
2123         }
2124         BOOST_UBLAS_INLINE
2125         iterator end () {
2126             return find (size_);
2127         }
2128
2129         // Reverse iterator
2130         typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
2131         typedef reverse_iterator_base<iterator> reverse_iterator;
2132
2133         BOOST_UBLAS_INLINE
2134         const_reverse_iterator rbegin () const {
2135             return const_reverse_iterator (end ());
2136         }
2137         BOOST_UBLAS_INLINE
2138         const_reverse_iterator crbegin () const {
2139             return rbegin ();
2140         }
2141         BOOST_UBLAS_INLINE
2142         const_reverse_iterator rend () const {
2143             return const_reverse_iterator (begin ());
2144         }
2145         BOOST_UBLAS_INLINE
2146         const_reverse_iterator crend () const {
2147             return rend ();
2148         }
2149         BOOST_UBLAS_INLINE
2150         reverse_iterator rbegin () {
2151             return reverse_iterator (end ());
2152         }
2153         BOOST_UBLAS_INLINE
2154         reverse_iterator rend () {
2155             return reverse_iterator (begin ());
2156         }
2157
2158          // Serialization
2159         template<class Archive>
2160         void serialize(Archive & ar, const unsigned int /* file_version */){
2161             serialization::collection_size_type s (size_);
2162             ar & serialization::make_nvp("size",s);
2163             if (Archive::is_loading::value) {
2164                 size_ = s;
2165             }
2166             // ISSUE: filled may be much less than capacity
2167             // ISSUE: index_data_ and value_data_ are undefined between filled and capacity (trouble with 'nan'-values)
2168             ar & serialization::make_nvp("capacity", capacity_);
2169             ar & serialization::make_nvp("filled", filled_);
2170             ar & serialization::make_nvp("sorted_filled", sorted_filled_);
2171             ar & serialization::make_nvp("sorted", sorted_);
2172             ar & serialization::make_nvp("index_data", index_data_);
2173             ar & serialization::make_nvp("value_data", value_data_);
2174             storage_invariants();
2175         }
2176
2177     private:
2178         void storage_invariants () const
2179         {
2180             BOOST_UBLAS_CHECK (capacity_ == index_data_.size (), internal_logic ());
2181             BOOST_UBLAS_CHECK (capacity_ == value_data_.size (), internal_logic ());
2182             BOOST_UBLAS_CHECK (filled_ <= capacity_, internal_logic ());
2183             BOOST_UBLAS_CHECK (sorted_filled_ <= filled_, internal_logic ());
2184             BOOST_UBLAS_CHECK (sorted_ == (sorted_filled_ == filled_), internal_logic ());
2185             BOOST_UBLAS_CHECK ((0 == filled_) || (zero_based(index_data_[filled_ - 1]) < size_), internal_logic ());
2186         }
2187
2188         size_type size_;
2189         size_type capacity_;
2190         mutable typename index_array_type::size_type filled_;
2191         mutable typename index_array_type::size_type sorted_filled_;
2192         mutable bool sorted_;
2193         mutable index_array_type index_data_;
2194         mutable value_array_type value_data_;
2195         static const value_type zero_;
2196
2197         BOOST_UBLAS_INLINE
2198         static size_type zero_based (size_type k_based_index) {
2199             return k_based_index - IB;
2200         }
2201         BOOST_UBLAS_INLINE
2202         static size_type k_based (size_type zero_based_index) {
2203             return zero_based_index + IB;
2204         }
2205
2206         friend class iterator;
2207         friend class const_iterator;
2208     };
2209
2210     template<class T, std::size_t IB, class IA, class TA>
2211     const typename coordinate_vector<T, IB, IA, TA>::value_type coordinate_vector<T, IB, IA, TA>::zero_ = value_type/*zero*/();
2212
2213 }}}
2214
2215 #endif