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     // Index map based sparse vector class
 264     template<class T, class A>
 265     class mapped_vector:
 266         public vector_container<mapped_vector<T, A> > {
 267
 268         typedef T &true_reference;
 269         typedef T *pointer;
 270         typedef const T *const_pointer;
 271         typedef mapped_vector<T, A> self_type;
 272     public:
 273 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
 274         using vector_container<self_type>::operator ();
 275 #endif
 276         typedef typename A::size_type size_type;
 277         typedef typename A::difference_type difference_type;
 278         typedef T value_type;
 279         typedef A array_type;
 280         typedef const value_type &const_reference;
 281 #ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
 282         typedef typename detail::map_traits<A,T>::reference reference;
 283 #else
 284         typedef sparse_vector_element<self_type> reference;
 285 #endif
 286         typedef const vector_reference<const self_type> const_closure_type;
 287         typedef vector_reference<self_type> closure_type;
 288         typedef self_type vector_temporary_type;
 289         typedef sparse_tag storage_category;
 290
 291         // Construction and destruction
 292         BOOST_UBLAS_INLINE
 293         mapped_vector ():
 294             vector_container<self_type> (),
 295             size_ (0), data_ () {}
 296         BOOST_UBLAS_INLINE
 297         mapped_vector (size_type size, size_type non_zeros = 0):
 298             vector_container<self_type> (),
 299             size_ (size), data_ () {
 300             detail::map_reserve (data(), restrict_capacity (non_zeros));
 301         }
 302         BOOST_UBLAS_INLINE
 303         mapped_vector (const mapped_vector &v):
 304             vector_container<self_type> (),
 305             size_ (v.size_), data_ (v.data_) {}
 306         template<class AE>
 307         BOOST_UBLAS_INLINE
 308         mapped_vector (const vector_expression<AE> &ae, size_type non_zeros = 0):
 309             vector_container<self_type> (),
 310             size_ (ae ().size ()), data_ () {
 311             detail::map_reserve (data(), restrict_capacity (non_zeros));
 312             vector_assign<scalar_assign> (*this, ae);
 313         }
 314
 315         // Accessors
 316         BOOST_UBLAS_INLINE
 317         size_type size () const {
 318             return size_;
 319         }
 320         BOOST_UBLAS_INLINE
 321         size_type nnz_capacity () const {
 322             return detail::map_capacity (data ());
 323         }
 324         BOOST_UBLAS_INLINE
 325         size_type nnz () const {
 326             return data (). size ();
 327         }
 328
 329         // Storage accessors
 330         BOOST_UBLAS_INLINE
 331         const array_type &data () const {
 332             return data_;
 333         }
 334         BOOST_UBLAS_INLINE
 335         array_type &data () {
 336             return data_;
 337         }
 338
 339         // Resizing
 340     private:
 341         BOOST_UBLAS_INLINE
 342         size_type restrict_capacity (size_type non_zeros) const {
 343             non_zeros = (std::min) (non_zeros, size_);
 344             return non_zeros;
 345         }
 346     public:
 347         BOOST_UBLAS_INLINE
 348         void resize (size_type size, bool preserve = true) {
 349             size_ = size;
 350             if (preserve) {
 351                 data ().erase (data ().lower_bound(size_), data ().end());
 352             }
 353             else {
 354                 data ().clear ();
 355             }
 356         }
 357
 358         // Reserving
 359         BOOST_UBLAS_INLINE
 360         void reserve (size_type non_zeros = 0, bool preserve = true) {
 361             detail::map_reserve (data (), restrict_capacity (non_zeros));
 362         }
 363
 364         // Element support
 365         BOOST_UBLAS_INLINE
 366         pointer find_element (size_type i) {
 367             return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i));
 368         }
 369         BOOST_UBLAS_INLINE
 370         const_pointer find_element (size_type i) const {
 371             const_subiterator_type it (data ().find (i));
 372             if (it == data ().end ())
 373                 return 0;
 374             BOOST_UBLAS_CHECK ((*it).first == i, internal_logic ());   // broken map
 375             return &(*it).second;
 376         }
 377
 378         // Element access
 379         BOOST_UBLAS_INLINE
 380         const_reference operator () (size_type i) const {
 381             BOOST_UBLAS_CHECK (i < size_, bad_index ());
 382             const_subiterator_type it (data ().find (i));
 383             if (it == data ().end ())
 384                 return zero_;
 385             BOOST_UBLAS_CHECK ((*it).first == i, internal_logic ());   // broken map
 386             return (*it).second;
 387         }
 388         BOOST_UBLAS_INLINE
 389         true_reference ref (size_type i) {
 390             BOOST_UBLAS_CHECK (i < size_, bad_index ());
 391             std::pair<subiterator_type, bool> ii (data ().insert (typename array_type::value_type (i, value_type/*zero*/())));
 392             BOOST_UBLAS_CHECK ((ii.first)->first == i, internal_logic ());   // broken map
 393             return (ii.first)->second;
 394         }
 395         BOOST_UBLAS_INLINE
 396         reference operator () (size_type i) {
 397 #ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
 398             return ref (i);
 399 #else
 400             BOOST_UBLAS_CHECK (i < size_, bad_index ());
 401             return reference (*this, i);
 402 #endif
 403         }
 404
 405         BOOST_UBLAS_INLINE
 406         const_reference operator [] (size_type i) const {
 407             return (*this) (i);
 408         }
 409         BOOST_UBLAS_INLINE
 410         reference operator [] (size_type i) {
 411             return (*this) (i);
 412         }
 413
 414         // Element assignment
 415         BOOST_UBLAS_INLINE
 416         true_reference insert_element (size_type i, const_reference t) {
 417             std::pair<subiterator_type, bool> ii = data ().insert (typename array_type::value_type (i, t));
 418             BOOST_UBLAS_CHECK (ii.second, bad_index ());        // duplicate element
 419             BOOST_UBLAS_CHECK ((ii.first)->first == i, internal_logic ());   // broken map
 420             if (!ii.second)     // existing element
 421                 (ii.first)->second = t;
 422             return (ii.first)->second;
 423         }
 424         BOOST_UBLAS_INLINE
 425         void erase_element (size_type i) {
 426             subiterator_type it = data ().find (i);
 427             if (it == data ().end ())
 428                 return;
 429             data ().erase (it);
 430         }
 431
 432         // Zeroing
 433         BOOST_UBLAS_INLINE
 434         void clear () {
 435             data ().clear ();
 436         }
 437
 438         // Assignment
 439         BOOST_UBLAS_INLINE
 440         mapped_vector &operator = (const mapped_vector &v) {
 441             if (this != &v) {
 442                 size_ = v.size_;
 443                 data () = v.data ();
 444             }
 445             return *this;
 446         }
 447         template<class C>          // Container assignment without temporary
 448         BOOST_UBLAS_INLINE
 449         mapped_vector &operator = (const vector_container<C> &v) {
 450             resize (v ().size (), false);
 451             assign (v);
 452             return *this;
 453         }
 454         BOOST_UBLAS_INLINE
 455         mapped_vector &assign_temporary (mapped_vector &v) {
 456             swap (v);
 457             return *this;
 458         }
 459         template<class AE>
 460         BOOST_UBLAS_INLINE
 461         mapped_vector &operator = (const vector_expression<AE> &ae) {
 462             self_type temporary (ae, detail::map_capacity (data()));
 463             return assign_temporary (temporary);
 464         }
 465         template<class AE>
 466         BOOST_UBLAS_INLINE
 467         mapped_vector &assign (const vector_expression<AE> &ae) {
 468             vector_assign<scalar_assign> (*this, ae);
 469             return *this;
 470         }
 471
 472         // Computed assignment
 473         template<class AE>
 474         BOOST_UBLAS_INLINE
 475         mapped_vector &operator += (const vector_expression<AE> &ae) {
 476             self_type temporary (*this + ae, detail::map_capacity (data()));
 477             return assign_temporary (temporary);
 478         }
 479         template<class C>          // Container assignment without temporary
 480         BOOST_UBLAS_INLINE
 481         mapped_vector &operator += (const vector_container<C> &v) {
 482             plus_assign (v);
 483             return *this;
 484         }
 485         template<class AE>
 486         BOOST_UBLAS_INLINE
 487         mapped_vector &plus_assign (const vector_expression<AE> &ae) {
 488             vector_assign<scalar_plus_assign> (*this, ae);
 489             return *this;
 490         }
 491         template<class AE>
 492         BOOST_UBLAS_INLINE
 493         mapped_vector &operator -= (const vector_expression<AE> &ae) {
 494             self_type temporary (*this - ae, detail::map_capacity (data()));
 495             return assign_temporary (temporary);
 496         }
 497         template<class C>          // Container assignment without temporary
 498         BOOST_UBLAS_INLINE
 499         mapped_vector &operator -= (const vector_container<C> &v) {
 500             minus_assign (v);
 501             return *this;
 502         }
 503         template<class AE>
 504         BOOST_UBLAS_INLINE
 505         mapped_vector &minus_assign (const vector_expression<AE> &ae) {
 506             vector_assign<scalar_minus_assign> (*this, ae);
 507             return *this;
 508         }
 509         template<class AT>
 510         BOOST_UBLAS_INLINE
 511         mapped_vector &operator *= (const AT &at) {
 512             vector_assign_scalar<scalar_multiplies_assign> (*this, at);
 513             return *this;
 514         }
 515         template<class AT>
 516         BOOST_UBLAS_INLINE
 517         mapped_vector &operator /= (const AT &at) {
 518             vector_assign_scalar<scalar_divides_assign> (*this, at);
 519             return *this;
 520         }
 521
 522         // Swapping
 523         BOOST_UBLAS_INLINE
 524         void swap (mapped_vector &v) {
 525             if (this != &v) {
 526                 std::swap (size_, v.size_);
 527                 data ().swap (v.data ());
 528             }
 529         }
 530         BOOST_UBLAS_INLINE
 531         friend void swap (mapped_vector &v1, mapped_vector &v2) {
 532             v1.swap (v2);
 533         }
 534
 535         // Iterator types
 536     private:
 537         // Use storage iterator
 538         typedef typename A::const_iterator const_subiterator_type;
 539         typedef typename A::iterator subiterator_type;
 540
 541         BOOST_UBLAS_INLINE
 542         true_reference at_element (size_type i) {
 543             BOOST_UBLAS_CHECK (i < size_, bad_index ());
 544             subiterator_type it (data ().find (i));
 545             BOOST_UBLAS_CHECK (it != data ().end(), bad_index ());
 546             BOOST_UBLAS_CHECK ((*it).first == i, internal_logic ());   // broken map
 547             return it->second;
 548         }
 549
 550     public:
 551         class const_iterator;
 552         class iterator;
 553
 554         // Element lookup
 555         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
 556         const_iterator find (size_type i) const {
 557             return const_iterator (*this, data ().lower_bound (i));
 558         }
 559         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
 560         iterator find (size_type i) {
 561             return iterator (*this, data ().lower_bound (i));
 562         }
 563
 564
 565         class const_iterator:
 566             public container_const_reference<mapped_vector>,
 567             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
 568                                                const_iterator, value_type> {
 569         public:
 570             typedef typename mapped_vector::value_type value_type;
 571             typedef typename mapped_vector::difference_type difference_type;
 572             typedef typename mapped_vector::const_reference reference;
 573             typedef const typename mapped_vector::pointer pointer;
 574
 575             // Construction and destruction
 576             BOOST_UBLAS_INLINE
 577             const_iterator ():
 578                 container_const_reference<self_type> (), it_ () {}
 579             BOOST_UBLAS_INLINE
 580             const_iterator (const self_type &v, const const_subiterator_type &it):
 581                 container_const_reference<self_type> (v), it_ (it) {}
 582             BOOST_UBLAS_INLINE
 583             const_iterator (const typename self_type::iterator &it):  // ISSUE self_type:: stops VC8 using std::iterator here
 584                 container_const_reference<self_type> (it ()), it_ (it.it_) {}
 585
 586             // Arithmetic
 587             BOOST_UBLAS_INLINE
 588             const_iterator &operator ++ () {
 589                 ++ it_;
 590                 return *this;
 591             }
 592             BOOST_UBLAS_INLINE
 593             const_iterator &operator -- () {
 594                 -- it_;
 595                 return *this;
 596             }
 597
 598             // Dereference
 599             BOOST_UBLAS_INLINE
 600             const_reference operator * () const {
 601                 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
 602                 return (*it_).second;
 603             }
 604
 605             // Index
 606             BOOST_UBLAS_INLINE
 607             size_type index () const {
 608                 BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
 609                 BOOST_UBLAS_CHECK ((*it_).first < (*this) ().size (), bad_index ());
 610                 return (*it_).first;
 611             }
 612
 613             // Assignment
 614             BOOST_UBLAS_INLINE
 615             const_iterator &operator = (const const_iterator &it) {
 616                 container_const_reference<self_type>::assign (&it ());
 617                 it_ = it.it_;
 618                 return *this;
 619             }
 620
 621             // Comparison
 622             BOOST_UBLAS_INLINE
 623             bool operator == (const const_iterator &it) const {
 624                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
 625                 return it_ == it.it_;
 626             }
 627
 628         private:
 629             const_subiterator_type it_;
 630         };
 631
 632         BOOST_UBLAS_INLINE
 633         const_iterator begin () const {
 634             return const_iterator (*this, data ().begin ());
 635         }
 636         BOOST_UBLAS_INLINE
 637         const_iterator end () const {
 638             return const_iterator (*this, data ().end ());
 639         }
 640
 641         class iterator:
 642             public container_reference<mapped_vector>,
 643             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
 644                                                iterator, value_type> {
 645         public:
 646             typedef typename mapped_vector::value_type value_type;
 647             typedef typename mapped_vector::difference_type difference_type;
 648             typedef typename mapped_vector::true_reference reference;
 649             typedef typename mapped_vector::pointer pointer;
 650
 651             // Construction and destruction
 652             BOOST_UBLAS_INLINE
 653             iterator ():
 654                 container_reference<self_type> (), it_ () {}
 655             BOOST_UBLAS_INLINE
 656             iterator (self_type &v, const subiterator_type &it):
 657                 container_reference<self_type> (v), it_ (it) {}
 658
 659             // Arithmetic
 660             BOOST_UBLAS_INLINE
 661             iterator &operator ++ () {
 662                 ++ it_;
 663                 return *this;
 664             }
 665             BOOST_UBLAS_INLINE
 666             iterator &operator -- () {
 667                 -- it_;
 668                 return *this;
 669             }
 670
 671             // Dereference
 672             BOOST_UBLAS_INLINE
 673             reference operator * () const {
 674                 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
 675                 return (*it_).second;
 676             }
 677
 678             // Index
 679             BOOST_UBLAS_INLINE
 680             size_type index () const {
 681                 BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
 682                 BOOST_UBLAS_CHECK ((*it_).first < (*this) ().size (), bad_index ());
 683                 return (*it_).first;
 684             }
 685
 686             // Assignment
 687             BOOST_UBLAS_INLINE
 688             iterator &operator = (const iterator &it) {
 689                 container_reference<self_type>::assign (&it ());
 690                 it_ = it.it_;
 691                 return *this;
 692             }
 693
 694             // Comparison
 695             BOOST_UBLAS_INLINE
 696             bool operator == (const iterator &it) const {
 697                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
 698                 return it_ == it.it_;
 699             }
 700
 701         private:
 702             subiterator_type it_;
 703
 704             friend class const_iterator;
 705         };
 706
 707         BOOST_UBLAS_INLINE
 708         iterator begin () {
 709             return iterator (*this, data ().begin ());
 710         }
 711         BOOST_UBLAS_INLINE
 712         iterator end () {
 713             return iterator (*this, data ().end ());
 714         }
 715
 716         // Reverse iterator
 717         typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
 718         typedef reverse_iterator_base<iterator> reverse_iterator;
 719
 720         BOOST_UBLAS_INLINE
 721         const_reverse_iterator rbegin () const {
 722             return const_reverse_iterator (end ());
 723         }
 724         BOOST_UBLAS_INLINE
 725         const_reverse_iterator rend () const {
 726             return const_reverse_iterator (begin ());
 727         }
 728         BOOST_UBLAS_INLINE
 729         reverse_iterator rbegin () {
 730             return reverse_iterator (end ());
 731         }
 732         BOOST_UBLAS_INLINE
 733         reverse_iterator rend () {
 734             return reverse_iterator (begin ());
 735         }
 736
 737          // Serialization
 738         template<class Archive>
 739         void serialize(Archive & ar, const unsigned int /* file_version */){
 740             serialization::collection_size_type s (size_);
 741             ar & serialization::make_nvp("size",s);
 742             if (Archive::is_loading::value) {
 743                 size_ = s;
 744             }
 745             ar & serialization::make_nvp("data", data_);
 746         }
 747
 748     private:
 749         size_type size_;
 750         array_type data_;
 751         static const value_type zero_;
 752     };
 753
 754     template<class T, class A>
 755     const typename mapped_vector<T, A>::value_type mapped_vector<T, A>::zero_ = value_type/*zero*/();
 756
 757
 758     // Compressed array based sparse vector class
 759     // Thanks to Kresimir Fresl for extending this to cover different index bases.
 760     template<class T, std::size_t IB, class IA, class TA>
 761     class compressed_vector:
 762         public vector_container<compressed_vector<T, IB, IA, TA> > {
 763
 764         typedef T &true_reference;
 765         typedef T *pointer;
 766         typedef const T *const_pointer;
 767         typedef compressed_vector<T, IB, IA, TA> self_type;
 768     public:
 769 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
 770         using vector_container<self_type>::operator ();
 771 #endif
 772         // ISSUE require type consistency check
 773         // is_convertable (IA::size_type, TA::size_type)
 774         typedef typename IA::value_type size_type;
 775         typedef typename IA::difference_type difference_type;
 776         typedef T value_type;
 777         typedef const T &const_reference;
 778 #ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
 779         typedef T &reference;
 780 #else
 781         typedef sparse_vector_element<self_type> reference;
 782 #endif
 783         typedef IA index_array_type;
 784         typedef TA value_array_type;
 785         typedef const vector_reference<const self_type> const_closure_type;
 786         typedef vector_reference<self_type> closure_type;
 787         typedef self_type vector_temporary_type;
 788         typedef sparse_tag storage_category;
 789
 790         // Construction and destruction
 791         BOOST_UBLAS_INLINE
 792         compressed_vector ():
 793             vector_container<self_type> (),
 794             size_ (0), capacity_ (restrict_capacity (0)), filled_ (0),
 795             index_data_ (capacity_), value_data_ (capacity_) {
 796             storage_invariants ();
 797         }
 798         explicit BOOST_UBLAS_INLINE
 799         compressed_vector (size_type size, size_type non_zeros = 0):
 800             vector_container<self_type> (),
 801             size_ (size), capacity_ (restrict_capacity (non_zeros)), filled_ (0),
 802             index_data_ (capacity_), value_data_ (capacity_) {
 803         storage_invariants ();
 804         }
 805         BOOST_UBLAS_INLINE
 806         compressed_vector (const compressed_vector &v):
 807             vector_container<self_type> (),
 808             size_ (v.size_), capacity_ (v.capacity_), filled_ (v.filled_),
 809             index_data_ (v.index_data_), value_data_ (v.value_data_) {
 810             storage_invariants ();
 811         }
 812         template<class AE>
 813         BOOST_UBLAS_INLINE
 814         compressed_vector (const vector_expression<AE> &ae, size_type non_zeros = 0):
 815             vector_container<self_type> (),
 816             size_ (ae ().size ()), capacity_ (restrict_capacity (non_zeros)), filled_ (0),
 817             index_data_ (capacity_), value_data_ (capacity_) {
 818             storage_invariants ();
 819             vector_assign<scalar_assign> (*this, ae);
 820         }
 821
 822         // Accessors
 823         BOOST_UBLAS_INLINE
 824         size_type size () const {
 825             return size_;
 826         }
 827         BOOST_UBLAS_INLINE
 828         size_type nnz_capacity () const {
 829             return capacity_;
 830         }
 831         BOOST_UBLAS_INLINE
 832         size_type nnz () const {
 833             return filled_;
 834         }
 835
 836         // Storage accessors
 837         BOOST_UBLAS_INLINE
 838         static size_type index_base () {
 839             return IB;
 840         }
 841         BOOST_UBLAS_INLINE
 842         typename index_array_type::size_type filled () const {
 843             return filled_;
 844         }
 845         BOOST_UBLAS_INLINE
 846         const index_array_type &index_data () const {
 847             return index_data_;
 848         }
 849         BOOST_UBLAS_INLINE
 850         const value_array_type &value_data () const {
 851             return value_data_;
 852         }
 853         BOOST_UBLAS_INLINE
 854         void set_filled (const typename index_array_type::size_type & filled) {
 855             filled_ = filled;
 856             storage_invariants ();
 857         }
 858         BOOST_UBLAS_INLINE
 859         index_array_type &index_data () {
 860             return index_data_;
 861         }
 862         BOOST_UBLAS_INLINE
 863         value_array_type &value_data () {
 864             return value_data_;
 865         }
 866
 867         // Resizing
 868     private:
 869         BOOST_UBLAS_INLINE
 870         size_type restrict_capacity (size_type non_zeros) const {
 871             non_zeros = (std::max) (non_zeros, size_type (1));
 872             non_zeros = (std::min) (non_zeros, size_);
 873             return non_zeros;
 874         }
 875     public:
 876         BOOST_UBLAS_INLINE
 877         void resize (size_type size, bool preserve = true) {
 878             size_ = size;
 879             capacity_ = restrict_capacity (capacity_);
 880             if (preserve) {
 881                 index_data_. resize (capacity_, size_type ());
 882                 value_data_. resize (capacity_, value_type ());
 883                 filled_ = (std::min) (capacity_, filled_);
 884                 while ((filled_ > 0) && (zero_based(index_data_[filled_ - 1]) >= size)) {
 885                     --filled_;
 886                 }
 887             }
 888             else {
 889                 index_data_. resize (capacity_);
 890                 value_data_. resize (capacity_);
 891                 filled_ = 0;
 892             }
 893             storage_invariants ();
 894         }
 895
 896         // Reserving
 897         BOOST_UBLAS_INLINE
 898         void reserve (size_type non_zeros, bool preserve = true) {
 899             capacity_ = restrict_capacity (non_zeros);
 900             if (preserve) {
 901                 index_data_. resize (capacity_, size_type ());
 902                 value_data_. resize (capacity_, value_type ());
 903                 filled_ = (std::min) (capacity_, filled_);
 904             }
 905             else {
 906                 index_data_. resize (capacity_);
 907                 value_data_. resize (capacity_);
 908                 filled_ = 0;
 909             }
 910             storage_invariants ();
 911         }
 912
 913         // Element support
 914         BOOST_UBLAS_INLINE
 915         pointer find_element (size_type i) {
 916             return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i));
 917         }
 918         BOOST_UBLAS_INLINE
 919         const_pointer find_element (size_type i) const {
 920             const_subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
 921             if (it == index_data_.begin () + filled_ || *it != k_based (i))
 922                 return 0;
 923             return &value_data_ [it - index_data_.begin ()];
 924         }
 925
 926         // Element access
 927         BOOST_UBLAS_INLINE
 928         const_reference operator () (size_type i) const {
 929             BOOST_UBLAS_CHECK (i < size_, bad_index ());
 930             const_subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
 931             if (it == index_data_.begin () + filled_ || *it != k_based (i))
 932                 return zero_;
 933             return value_data_ [it - index_data_.begin ()];
 934         }
 935         BOOST_UBLAS_INLINE
 936         true_reference ref (size_type i) {
 937             BOOST_UBLAS_CHECK (i < size_, bad_index ());
 938             subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
 939             if (it == index_data_.begin () + filled_ || *it != k_based (i))
 940                 return insert_element (i, value_type/*zero*/());
 941             else
 942                 return value_data_ [it - index_data_.begin ()];
 943         }
 944         BOOST_UBLAS_INLINE
 945         reference operator () (size_type i) {
 946 #ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
 947             return ref (i) ;
 948 #else
 949             BOOST_UBLAS_CHECK (i < size_, bad_index ());
 950             return reference (*this, i);
 951 #endif
 952         }
 953
 954         BOOST_UBLAS_INLINE
 955         const_reference operator [] (size_type i) const {
 956             return (*this) (i);
 957         }
 958         BOOST_UBLAS_INLINE
 959         reference operator [] (size_type i) {
 960             return (*this) (i);
 961         }
 962
 963         // Element assignment
 964         BOOST_UBLAS_INLINE
 965         true_reference insert_element (size_type i, const_reference t) {
 966             BOOST_UBLAS_CHECK (!find_element (i), bad_index ());        // duplicate element
 967             if (filled_ >= capacity_)
 968                 reserve (2 * capacity_, true);
 969             subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
 970             // ISSUE max_capacity limit due to difference_type
 971             typename std::iterator_traits<subiterator_type>::difference_type n = it - index_data_.begin ();
 972             BOOST_UBLAS_CHECK (filled_ == 0 || filled_ == typename index_array_type::size_type (n) || *it != k_based (i), internal_logic ());   // duplicate found by lower_bound
 973             ++ filled_;
 974             it = index_data_.begin () + n;
 975             std::copy_backward (it, index_data_.begin () + filled_ - 1, index_data_.begin () + filled_);
 976             *it = k_based (i);
 977             typename value_array_type::iterator itt (value_data_.begin () + n);
 978             std::copy_backward (itt, value_data_.begin () + filled_ - 1, value_data_.begin () + filled_);
 979             *itt = t;
 980             storage_invariants ();
 981             return *itt;
 982         }
 983         BOOST_UBLAS_INLINE
 984         void erase_element (size_type i) {
 985             subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
 986             typename std::iterator_traits<subiterator_type>::difference_type  n = it - index_data_.begin ();
 987             if (filled_ > typename index_array_type::size_type (n) && *it == k_based (i)) {
 988                 std::copy (it + 1, index_data_.begin () + filled_, it);
 989                 typename value_array_type::iterator itt (value_data_.begin () + n);
 990                 std::copy (itt + 1, value_data_.begin () + filled_, itt);
 991                 -- filled_;
 992             }
 993             storage_invariants ();
 994         }
 995
 996         // Zeroing
 997         BOOST_UBLAS_INLINE
 998         void clear () {
 999             filled_ = 0;
1000             storage_invariants ();
1001         }
1002
1003         // Assignment
1004         BOOST_UBLAS_INLINE
1005         compressed_vector &operator = (const compressed_vector &v) {
1006             if (this != &v) {
1007                 size_ = v.size_;
1008                 capacity_ = v.capacity_;
1009                 filled_ = v.filled_;
1010                 index_data_ = v.index_data_;
1011                 value_data_ = v.value_data_;
1012             }
1013             storage_invariants ();
1014             return *this;
1015         }
1016         template<class C>          // Container assignment without temporary
1017         BOOST_UBLAS_INLINE
1018         compressed_vector &operator = (const vector_container<C> &v) {
1019             resize (v ().size (), false);
1020             assign (v);
1021             return *this;
1022         }
1023         BOOST_UBLAS_INLINE
1024         compressed_vector &assign_temporary (compressed_vector &v) {
1025             swap (v);
1026             return *this;
1027         }
1028         template<class AE>
1029         BOOST_UBLAS_INLINE
1030         compressed_vector &operator = (const vector_expression<AE> &ae) {
1031             self_type temporary (ae, capacity_);
1032             return assign_temporary (temporary);
1033         }
1034         template<class AE>
1035         BOOST_UBLAS_INLINE
1036         compressed_vector &assign (const vector_expression<AE> &ae) {
1037             vector_assign<scalar_assign> (*this, ae);
1038             return *this;
1039         }
1040
1041         // Computed assignment
1042         template<class AE>
1043         BOOST_UBLAS_INLINE
1044         compressed_vector &operator += (const vector_expression<AE> &ae) {
1045             self_type temporary (*this + ae, capacity_);
1046             return assign_temporary (temporary);
1047         }
1048         template<class C>          // Container assignment without temporary
1049         BOOST_UBLAS_INLINE
1050         compressed_vector &operator += (const vector_container<C> &v) {
1051             plus_assign (v);
1052             return *this;
1053         }
1054         template<class AE>
1055         BOOST_UBLAS_INLINE
1056         compressed_vector &plus_assign (const vector_expression<AE> &ae) {
1057             vector_assign<scalar_plus_assign> (*this, ae);
1058             return *this;
1059         }
1060         template<class AE>
1061         BOOST_UBLAS_INLINE
1062         compressed_vector &operator -= (const vector_expression<AE> &ae) {
1063             self_type temporary (*this - ae, capacity_);
1064             return assign_temporary (temporary);
1065         }
1066         template<class C>          // Container assignment without temporary
1067         BOOST_UBLAS_INLINE
1068         compressed_vector &operator -= (const vector_container<C> &v) {
1069             minus_assign (v);
1070             return *this;
1071         }
1072         template<class AE>
1073         BOOST_UBLAS_INLINE
1074         compressed_vector &minus_assign (const vector_expression<AE> &ae) {
1075             vector_assign<scalar_minus_assign> (*this, ae);
1076             return *this;
1077         }
1078         template<class AT>
1079         BOOST_UBLAS_INLINE
1080         compressed_vector &operator *= (const AT &at) {
1081             vector_assign_scalar<scalar_multiplies_assign> (*this, at);
1082             return *this;
1083         }
1084         template<class AT>
1085         BOOST_UBLAS_INLINE
1086         compressed_vector &operator /= (const AT &at) {
1087             vector_assign_scalar<scalar_divides_assign> (*this, at);
1088             return *this;
1089         }
1090
1091         // Swapping
1092         BOOST_UBLAS_INLINE
1093         void swap (compressed_vector &v) {
1094             if (this != &v) {
1095                 std::swap (size_, v.size_);
1096                 std::swap (capacity_, v.capacity_);
1097                 std::swap (filled_, v.filled_);
1098                 index_data_.swap (v.index_data_);
1099                 value_data_.swap (v.value_data_);
1100             }
1101             storage_invariants ();
1102         }
1103         BOOST_UBLAS_INLINE
1104         friend void swap (compressed_vector &v1, compressed_vector &v2) {
1105             v1.swap (v2);
1106         }
1107
1108         // Back element insertion and erasure
1109         BOOST_UBLAS_INLINE
1110         void push_back (size_type i, const_reference t) {
1111             BOOST_UBLAS_CHECK (filled_ == 0 || index_data_ [filled_ - 1] < k_based (i), external_logic ());
1112             if (filled_ >= capacity_)
1113                 reserve (2 * capacity_, true);
1114             BOOST_UBLAS_CHECK (filled_ < capacity_, internal_logic ());
1115             index_data_ [filled_] = k_based (i);
1116             value_data_ [filled_] = t;
1117             ++ filled_;
1118             storage_invariants ();
1119         }
1120         BOOST_UBLAS_INLINE
1121         void pop_back () {
1122             BOOST_UBLAS_CHECK (filled_ > 0, external_logic ());
1123             -- filled_;
1124             storage_invariants ();
1125         }
1126
1127         // Iterator types
1128     private:
1129         // Use index array iterator
1130         typedef typename IA::const_iterator const_subiterator_type;
1131         typedef typename IA::iterator subiterator_type;
1132
1133         BOOST_UBLAS_INLINE
1134         true_reference at_element (size_type i) {
1135             BOOST_UBLAS_CHECK (i < size_, bad_index ());
1136             subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1137             BOOST_UBLAS_CHECK (it != index_data_.begin () + filled_ && *it == k_based (i), bad_index ());
1138             return value_data_ [it - index_data_.begin ()];
1139         }
1140
1141     public:
1142         class const_iterator;
1143         class iterator;
1144
1145         // Element lookup
1146         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1147         const_iterator find (size_type i) const {
1148             return const_iterator (*this, detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1149         }
1150         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1151         iterator find (size_type i) {
1152             return iterator (*this, detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1153         }
1154
1155
1156         class const_iterator:
1157             public container_const_reference<compressed_vector>,
1158             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
1159                                                const_iterator, value_type> {
1160         public:
1161             typedef typename compressed_vector::value_type value_type;
1162             typedef typename compressed_vector::difference_type difference_type;
1163             typedef typename compressed_vector::const_reference reference;
1164             typedef const typename compressed_vector::pointer pointer;
1165
1166             // Construction and destruction
1167             BOOST_UBLAS_INLINE
1168             const_iterator ():
1169                 container_const_reference<self_type> (), it_ () {}
1170             BOOST_UBLAS_INLINE
1171             const_iterator (const self_type &v, const const_subiterator_type &it):
1172                 container_const_reference<self_type> (v), it_ (it) {}
1173             BOOST_UBLAS_INLINE
1174             const_iterator (const typename self_type::iterator &it):  // ISSUE self_type:: stops VC8 using std::iterator here
1175                 container_const_reference<self_type> (it ()), it_ (it.it_) {}
1176
1177             // Arithmetic
1178             BOOST_UBLAS_INLINE
1179             const_iterator &operator ++ () {
1180                 ++ it_;
1181                 return *this;
1182             }
1183             BOOST_UBLAS_INLINE
1184             const_iterator &operator -- () {
1185                 -- it_;
1186                 return *this;
1187             }
1188
1189             // Dereference
1190             BOOST_UBLAS_INLINE
1191             const_reference operator * () const {
1192                 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
1193                 return (*this) ().value_data_ [it_ - (*this) ().index_data_.begin ()];
1194             }
1195
1196             // Index
1197             BOOST_UBLAS_INLINE
1198             size_type index () const {
1199                 BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
1200                 BOOST_UBLAS_CHECK ((*this) ().zero_based (*it_) < (*this) ().size (), bad_index ());
1201                 return (*this) ().zero_based (*it_);
1202             }
1203
1204             // Assignment
1205             BOOST_UBLAS_INLINE
1206             const_iterator &operator = (const const_iterator &it) {
1207                 container_const_reference<self_type>::assign (&it ());
1208                 it_ = it.it_;
1209                 return *this;
1210             }
1211
1212             // Comparison
1213             BOOST_UBLAS_INLINE
1214             bool operator == (const const_iterator &it) const {
1215                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1216                 return it_ == it.it_;
1217             }
1218
1219         private:
1220             const_subiterator_type it_;
1221         };
1222
1223         BOOST_UBLAS_INLINE
1224         const_iterator begin () const {
1225             return find (0);
1226         }
1227         BOOST_UBLAS_INLINE
1228         const_iterator end () const {
1229             return find (size_);
1230         }
1231
1232         class iterator:
1233             public container_reference<compressed_vector>,
1234             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
1235                                                iterator, value_type> {
1236         public:
1237             typedef typename compressed_vector::value_type value_type;
1238             typedef typename compressed_vector::difference_type difference_type;
1239             typedef typename compressed_vector::true_reference reference;
1240             typedef typename compressed_vector::pointer pointer;
1241
1242             // Construction and destruction
1243             BOOST_UBLAS_INLINE
1244             iterator ():
1245                 container_reference<self_type> (), it_ () {}
1246             BOOST_UBLAS_INLINE
1247             iterator (self_type &v, const subiterator_type &it):
1248                 container_reference<self_type> (v), it_ (it) {}
1249
1250             // Arithmetic
1251             BOOST_UBLAS_INLINE
1252             iterator &operator ++ () {
1253                 ++ it_;
1254                 return *this;
1255             }
1256             BOOST_UBLAS_INLINE
1257             iterator &operator -- () {
1258                 -- it_;
1259                 return *this;
1260             }
1261
1262             // Dereference
1263             BOOST_UBLAS_INLINE
1264             reference operator * () const {
1265                 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
1266                 return (*this) ().value_data_ [it_ - (*this) ().index_data_.begin ()];
1267             }
1268
1269             // Index
1270             BOOST_UBLAS_INLINE
1271             size_type index () const {
1272                 BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
1273                 BOOST_UBLAS_CHECK ((*this) ().zero_based (*it_) < (*this) ().size (), bad_index ());
1274                 return (*this) ().zero_based (*it_);
1275             }
1276
1277             // Assignment
1278             BOOST_UBLAS_INLINE
1279             iterator &operator = (const iterator &it) {
1280                 container_reference<self_type>::assign (&it ());
1281                 it_ = it.it_;
1282                 return *this;
1283             }
1284
1285             // Comparison
1286             BOOST_UBLAS_INLINE
1287             bool operator == (const iterator &it) const {
1288                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1289                 return it_ == it.it_;
1290             }
1291
1292         private:
1293             subiterator_type it_;
1294
1295             friend class const_iterator;
1296         };
1297
1298         BOOST_UBLAS_INLINE
1299         iterator begin () {
1300             return find (0);
1301         }
1302         BOOST_UBLAS_INLINE
1303         iterator end () {
1304             return find (size_);
1305         }
1306
1307         // Reverse iterator
1308         typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
1309         typedef reverse_iterator_base<iterator> reverse_iterator;
1310
1311         BOOST_UBLAS_INLINE
1312         const_reverse_iterator rbegin () const {
1313             return const_reverse_iterator (end ());
1314         }
1315         BOOST_UBLAS_INLINE
1316         const_reverse_iterator rend () const {
1317             return const_reverse_iterator (begin ());
1318         }
1319         BOOST_UBLAS_INLINE
1320         reverse_iterator rbegin () {
1321             return reverse_iterator (end ());
1322         }
1323         BOOST_UBLAS_INLINE
1324         reverse_iterator rend () {
1325             return reverse_iterator (begin ());
1326         }
1327
1328          // Serialization
1329         template<class Archive>
1330         void serialize(Archive & ar, const unsigned int /* file_version */){
1331             serialization::collection_size_type s (size_);
1332             ar & serialization::make_nvp("size",s);
1333             if (Archive::is_loading::value) {
1334                 size_ = s;
1335             }
1336             // ISSUE: filled may be much less than capacity
1337             // ISSUE: index_data_ and value_data_ are undefined between filled and capacity (trouble with 'nan'-values)
1338             ar & serialization::make_nvp("capacity", capacity_);
1339             ar & serialization::make_nvp("filled", filled_);
1340             ar & serialization::make_nvp("index_data", index_data_);
1341             ar & serialization::make_nvp("value_data", value_data_);
1342             storage_invariants();
1343         }
1344
1345     private:
1346         void storage_invariants () const
1347         {
1348             BOOST_UBLAS_CHECK (capacity_ == index_data_.size (), internal_logic ());
1349             BOOST_UBLAS_CHECK (capacity_ == value_data_.size (), internal_logic ());
1350             BOOST_UBLAS_CHECK (filled_ <= capacity_, internal_logic ());
1351             BOOST_UBLAS_CHECK ((0 == filled_) || (zero_based(index_data_[filled_ - 1]) < size_), internal_logic ());
1352         }
1353
1354         size_type size_;
1355         typename index_array_type::size_type capacity_;
1356         typename index_array_type::size_type filled_;
1357         index_array_type index_data_;
1358         value_array_type value_data_;
1359         static const value_type zero_;
1360
1361         BOOST_UBLAS_INLINE
1362         static size_type zero_based (size_type k_based_index) {
1363             return k_based_index - IB;
1364         }
1365         BOOST_UBLAS_INLINE
1366         static size_type k_based (size_type zero_based_index) {
1367             return zero_based_index + IB;
1368         }
1369
1370         friend class iterator;
1371         friend class const_iterator;
1372     };
1373
1374     template<class T, std::size_t IB, class IA, class TA>
1375     const typename compressed_vector<T, IB, IA, TA>::value_type compressed_vector<T, IB, IA, TA>::zero_ = value_type/*zero*/();
1376
1377
1378     // Coordimate array based sparse vector class
1379     // Thanks to Kresimir Fresl for extending this to cover different index bases.
1380     template<class T, std::size_t IB, class IA, class TA>
1381     class coordinate_vector:
1382         public vector_container<coordinate_vector<T, IB, IA, TA> > {
1383
1384         typedef T &true_reference;
1385         typedef T *pointer;
1386         typedef const T *const_pointer;
1387         typedef coordinate_vector<T, IB, IA, TA> self_type;
1388     public:
1389 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
1390         using vector_container<self_type>::operator ();
1391 #endif
1392         // ISSUE require type consistency check
1393         // is_convertable (IA::size_type, TA::size_type)
1394         typedef typename IA::value_type size_type;
1395         typedef typename IA::difference_type difference_type;
1396         typedef T value_type;
1397         typedef const T &const_reference;
1398 #ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
1399         typedef T &reference;
1400 #else
1401         typedef sparse_vector_element<self_type> reference;
1402 #endif
1403         typedef IA index_array_type;
1404         typedef TA value_array_type;
1405         typedef const vector_reference<const self_type> const_closure_type;
1406         typedef vector_reference<self_type> closure_type;
1407         typedef self_type vector_temporary_type;
1408         typedef sparse_tag storage_category;
1409
1410         // Construction and destruction
1411         BOOST_UBLAS_INLINE
1412         coordinate_vector ():
1413             vector_container<self_type> (),
1414             size_ (0), capacity_ (restrict_capacity (0)),
1415             filled_ (0), sorted_filled_ (filled_), sorted_ (true),
1416             index_data_ (capacity_), value_data_ (capacity_) {
1417             storage_invariants ();
1418         }
1419         explicit BOOST_UBLAS_INLINE
1420         coordinate_vector (size_type size, size_type non_zeros = 0):
1421             vector_container<self_type> (),
1422             size_ (size), capacity_ (restrict_capacity (non_zeros)),
1423             filled_ (0), sorted_filled_ (filled_), sorted_ (true),
1424             index_data_ (capacity_), value_data_ (capacity_) {
1425             storage_invariants ();
1426         }
1427         BOOST_UBLAS_INLINE
1428         coordinate_vector (const coordinate_vector &v):
1429             vector_container<self_type> (),
1430             size_ (v.size_), capacity_ (v.capacity_),
1431             filled_ (v.filled_), sorted_filled_ (v.sorted_filled_), sorted_ (v.sorted_),
1432             index_data_ (v.index_data_), value_data_ (v.value_data_) {
1433             storage_invariants ();
1434         }
1435         template<class AE>
1436         BOOST_UBLAS_INLINE
1437         coordinate_vector (const vector_expression<AE> &ae, size_type non_zeros = 0):
1438             vector_container<self_type> (),
1439             size_ (ae ().size ()), capacity_ (restrict_capacity (non_zeros)),
1440             filled_ (0), sorted_filled_ (filled_), sorted_ (true),
1441             index_data_ (capacity_), value_data_ (capacity_) {
1442             storage_invariants ();
1443             vector_assign<scalar_assign> (*this, ae);
1444         }
1445
1446         // Accessors
1447         BOOST_UBLAS_INLINE
1448         size_type size () const {
1449             return size_;
1450         }
1451         BOOST_UBLAS_INLINE
1452         size_type nnz_capacity () const {
1453             return capacity_;
1454         }
1455         BOOST_UBLAS_INLINE
1456         size_type nnz () const {
1457             return filled_;
1458         }
1459
1460         // Storage accessors
1461         BOOST_UBLAS_INLINE
1462         static size_type index_base () {
1463             return IB;
1464         }
1465         BOOST_UBLAS_INLINE
1466         typename index_array_type::size_type filled () const {
1467             return filled_;
1468         }
1469         BOOST_UBLAS_INLINE
1470         const index_array_type &index_data () const {
1471             return index_data_;
1472         }
1473         BOOST_UBLAS_INLINE
1474         const value_array_type &value_data () const {
1475             return value_data_;
1476         }
1477         BOOST_UBLAS_INLINE
1478         void set_filled (const typename index_array_type::size_type &sorted, const typename index_array_type::size_type &filled) {
1479             sorted_filled_ = sorted;
1480             filled_ = filled;
1481             storage_invariants ();
1482             return filled_;
1483         }
1484         BOOST_UBLAS_INLINE
1485         index_array_type &index_data () {
1486             return index_data_;
1487         }
1488         BOOST_UBLAS_INLINE
1489         value_array_type &value_data () {
1490             return value_data_;
1491         }
1492
1493         // Resizing
1494     private:
1495         BOOST_UBLAS_INLINE
1496         size_type restrict_capacity (size_type non_zeros) const {
1497              // minimum non_zeros
1498              non_zeros = (std::max) (non_zeros, size_type (1));
1499              // ISSUE no maximum as coordinate may contain inserted duplicates
1500              return non_zeros;
1501         }
1502     public:
1503         BOOST_UBLAS_INLINE
1504         void resize (size_type size, bool preserve = true) {
1505             if (preserve)
1506                 sort ();    // remove duplicate elements.
1507             size_ = size;
1508             capacity_ = restrict_capacity (capacity_);
1509             if (preserve) {
1510                 index_data_. resize (capacity_, size_type ());
1511                 value_data_. resize (capacity_, value_type ());
1512                 filled_ = (std::min) (capacity_, filled_);
1513                 while ((filled_ > 0) && (zero_based(index_data_[filled_ - 1]) >= size)) {
1514                     --filled_;
1515                 }
1516             }
1517             else {
1518                 index_data_. resize (capacity_);
1519                 value_data_. resize (capacity_);
1520                 filled_ = 0;
1521             }
1522             sorted_filled_ = filled_;
1523             storage_invariants ();
1524         }
1525         // Reserving
1526         BOOST_UBLAS_INLINE
1527         void reserve (size_type non_zeros, bool preserve = true) {
1528             if (preserve)
1529                 sort ();    // remove duplicate elements.
1530             capacity_ = restrict_capacity (non_zeros);
1531             if (preserve) {
1532                 index_data_. resize (capacity_, size_type ());
1533                 value_data_. resize (capacity_, value_type ());
1534                 filled_ = (std::min) (capacity_, filled_);
1535                 }
1536             else {
1537                 index_data_. resize (capacity_);
1538                 value_data_. resize (capacity_);
1539                 filled_ = 0;
1540             }
1541             sorted_filled_ = filled_;
1542             storage_invariants ();
1543         }
1544
1545         // Element support
1546         BOOST_UBLAS_INLINE
1547         pointer find_element (size_type i) {
1548             return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i));
1549         }
1550         BOOST_UBLAS_INLINE
1551         const_pointer find_element (size_type i) const {
1552             sort ();
1553             const_subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1554             if (it == index_data_.begin () + filled_ || *it != k_based (i))
1555                 return 0;
1556             return &value_data_ [it - index_data_.begin ()];
1557         }
1558
1559         // Element access
1560         BOOST_UBLAS_INLINE
1561         const_reference operator () (size_type i) const {
1562             BOOST_UBLAS_CHECK (i < size_, bad_index ());
1563             sort ();
1564             const_subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1565             if (it == index_data_.begin () + filled_ || *it != k_based (i))
1566                 return zero_;
1567             return value_data_ [it - index_data_.begin ()];
1568         }
1569         BOOST_UBLAS_INLINE
1570         true_reference ref (size_type i) {
1571             BOOST_UBLAS_CHECK (i < size_, bad_index ());
1572             sort ();
1573             subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1574             if (it == index_data_.begin () + filled_ || *it != k_based (i))
1575                 return insert_element (i, value_type/*zero*/());
1576             else
1577                 return value_data_ [it - index_data_.begin ()];
1578         }
1579         BOOST_UBLAS_INLINE
1580         reference operator () (size_type i) {
1581 #ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
1582             return ref (i);
1583 #else
1584             BOOST_UBLAS_CHECK (i < size_, bad_index ());
1585             return reference (*this, i);
1586 #endif
1587         }
1588
1589         BOOST_UBLAS_INLINE
1590         const_reference operator [] (size_type i) const {
1591             return (*this) (i);
1592         }
1593         BOOST_UBLAS_INLINE
1594         reference operator [] (size_type i) {
1595             return (*this) (i);
1596         }
1597
1598         // Element assignment
1599         BOOST_UBLAS_INLINE
1600         void append_element (size_type i, const_reference t) {
1601             if (filled_ >= capacity_)
1602                 reserve (2 * filled_, true);
1603             BOOST_UBLAS_CHECK (filled_ < capacity_, internal_logic ());
1604             index_data_ [filled_] = k_based (i);
1605             value_data_ [filled_] = t;
1606             ++ filled_;
1607             sorted_ = false;
1608             storage_invariants ();
1609         }
1610         BOOST_UBLAS_INLINE
1611         true_reference insert_element (size_type i, const_reference t) {
1612             BOOST_UBLAS_CHECK (!find_element (i), bad_index ());        // duplicate element
1613             append_element (i, t);
1614             return value_data_ [filled_ - 1];
1615         }
1616         BOOST_UBLAS_INLINE
1617         void erase_element (size_type i) {
1618             sort ();
1619             subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1620             typename std::iterator_traits<subiterator_type>::difference_type n = it - index_data_.begin ();
1621             if (filled_ > typename index_array_type::size_type (n) && *it == k_based (i)) {
1622                 std::copy (it + 1, index_data_.begin () + filled_, it);
1623                 typename value_array_type::iterator itt (value_data_.begin () + n);
1624                 std::copy (itt + 1, value_data_.begin () + filled_, itt);
1625                 -- filled_;
1626                 sorted_filled_ = filled_;
1627             }
1628             storage_invariants ();
1629         }
1630
1631         // Zeroing
1632         BOOST_UBLAS_INLINE
1633         void clear () {
1634             filled_ = 0;
1635             sorted_filled_ = filled_;
1636             sorted_ = true;
1637             storage_invariants ();
1638         }
1639
1640         // Assignment
1641         BOOST_UBLAS_INLINE
1642         coordinate_vector &operator = (const coordinate_vector &v) {
1643             if (this != &v) {
1644                 size_ = v.size_;
1645                 capacity_ = v.capacity_;
1646                 filled_ = v.filled_;
1647                 sorted_filled_ = v.sorted_filled_;
1648                 sorted_ = v.sorted_;
1649                 index_data_ = v.index_data_;
1650                 value_data_ = v.value_data_;
1651             }
1652             storage_invariants ();
1653             return *this;
1654         }
1655         template<class C>          // Container assignment without temporary
1656         BOOST_UBLAS_INLINE
1657         coordinate_vector &operator = (const vector_container<C> &v) {
1658             resize (v ().size (), false);
1659             assign (v);
1660             return *this;
1661         }
1662         BOOST_UBLAS_INLINE
1663         coordinate_vector &assign_temporary (coordinate_vector &v) {
1664             swap (v);
1665             return *this;
1666         }
1667         template<class AE>
1668         BOOST_UBLAS_INLINE
1669         coordinate_vector &operator = (const vector_expression<AE> &ae) {
1670             self_type temporary (ae, capacity_);
1671             return assign_temporary (temporary);
1672         }
1673         template<class AE>
1674         BOOST_UBLAS_INLINE
1675         coordinate_vector &assign (const vector_expression<AE> &ae) {
1676             vector_assign<scalar_assign> (*this, ae);
1677             return *this;
1678         }
1679
1680         // Computed assignment
1681         template<class AE>
1682         BOOST_UBLAS_INLINE
1683         coordinate_vector &operator += (const vector_expression<AE> &ae) {
1684             self_type temporary (*this + ae, capacity_);
1685             return assign_temporary (temporary);
1686         }
1687         template<class C>          // Container assignment without temporary
1688         BOOST_UBLAS_INLINE
1689         coordinate_vector &operator += (const vector_container<C> &v) {
1690             plus_assign (v);
1691             return *this;
1692         }
1693         template<class AE>
1694         BOOST_UBLAS_INLINE
1695         coordinate_vector &plus_assign (const vector_expression<AE> &ae) {
1696             vector_assign<scalar_plus_assign> (*this, ae);
1697             return *this;
1698         }
1699         template<class AE>
1700         BOOST_UBLAS_INLINE
1701         coordinate_vector &operator -= (const vector_expression<AE> &ae) {
1702             self_type temporary (*this - ae, capacity_);
1703             return assign_temporary (temporary);
1704         }
1705         template<class C>          // Container assignment without temporary
1706         BOOST_UBLAS_INLINE
1707         coordinate_vector &operator -= (const vector_container<C> &v) {
1708             minus_assign (v);
1709             return *this;
1710         }
1711         template<class AE>
1712         BOOST_UBLAS_INLINE
1713         coordinate_vector &minus_assign (const vector_expression<AE> &ae) {
1714             vector_assign<scalar_minus_assign> (*this, ae);
1715             return *this;
1716         }
1717         template<class AT>
1718         BOOST_UBLAS_INLINE
1719         coordinate_vector &operator *= (const AT &at) {
1720             vector_assign_scalar<scalar_multiplies_assign> (*this, at);
1721             return *this;
1722         }
1723         template<class AT>
1724         BOOST_UBLAS_INLINE
1725         coordinate_vector &operator /= (const AT &at) {
1726             vector_assign_scalar<scalar_divides_assign> (*this, at);
1727             return *this;
1728         }
1729
1730         // Swapping
1731         BOOST_UBLAS_INLINE
1732         void swap (coordinate_vector &v) {
1733             if (this != &v) {
1734                 std::swap (size_, v.size_);
1735                 std::swap (capacity_, v.capacity_);
1736                 std::swap (filled_, v.filled_);
1737                 std::swap (sorted_filled_, v.sorted_filled_);
1738                 std::swap (sorted_, v.sorted_);
1739                 index_data_.swap (v.index_data_);
1740                 value_data_.swap (v.value_data_);
1741             }
1742             storage_invariants ();
1743         }
1744         BOOST_UBLAS_INLINE
1745         friend void swap (coordinate_vector &v1, coordinate_vector &v2) {
1746             v1.swap (v2);
1747         }
1748
1749         // Sorting and summation of duplicates
1750         BOOST_UBLAS_INLINE
1751         void sort () const {
1752             if (! sorted_ && filled_ > 0) {
1753                 typedef index_pair_array<index_array_type, value_array_type> array_pair;
1754                 array_pair ipa (filled_, index_data_, value_data_);
1755                 const typename array_pair::iterator iunsorted = ipa.begin () + sorted_filled_;
1756                 // sort new elements and merge
1757                 std::sort (iunsorted, ipa.end ());
1758                 std::inplace_merge (ipa.begin (), iunsorted, ipa.end ());
1759
1760                 // sum duplicates with += and remove
1761                 size_type filled = 0;
1762                 for (size_type i = 1; i < filled_; ++ i) {
1763                     if (index_data_ [filled] != index_data_ [i]) {
1764                         ++ filled;
1765                         if (filled != i) {
1766                             index_data_ [filled] = index_data_ [i];
1767                             value_data_ [filled] = value_data_ [i];
1768                         }
1769                     } else {
1770                         value_data_ [filled] += value_data_ [i];
1771                     }
1772                 }
1773                 filled_ = filled + 1;
1774                 sorted_filled_ = filled_;
1775                 sorted_ = true;
1776                 storage_invariants ();
1777             }
1778         }
1779
1780         // Back element insertion and erasure
1781         BOOST_UBLAS_INLINE
1782         void push_back (size_type i, const_reference t) {
1783             // must maintain sort order
1784             BOOST_UBLAS_CHECK (sorted_ && (filled_ == 0 || index_data_ [filled_ - 1] < k_based (i)), external_logic ());
1785             if (filled_ >= capacity_)
1786                 reserve (2 * filled_, true);
1787             BOOST_UBLAS_CHECK (filled_ < capacity_, internal_logic ());
1788             index_data_ [filled_] = k_based (i);
1789             value_data_ [filled_] = t;
1790             ++ filled_;
1791             sorted_filled_ = filled_;
1792             storage_invariants ();
1793         }
1794         BOOST_UBLAS_INLINE
1795         void pop_back () {
1796             // ISSUE invariants could be simpilfied if sorted required as precondition
1797             BOOST_UBLAS_CHECK (filled_ > 0, external_logic ());
1798             -- filled_;
1799             sorted_filled_ = (std::min) (sorted_filled_, filled_);
1800             sorted_ = sorted_filled_ = filled_;
1801             storage_invariants ();
1802         }
1803
1804         // Iterator types
1805     private:
1806         // Use index array iterator
1807         typedef typename IA::const_iterator const_subiterator_type;
1808         typedef typename IA::iterator subiterator_type;
1809
1810         BOOST_UBLAS_INLINE
1811         true_reference at_element (size_type i) {
1812             BOOST_UBLAS_CHECK (i < size_, bad_index ());
1813             sort ();
1814             subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1815             BOOST_UBLAS_CHECK (it != index_data_.begin () + filled_ && *it == k_based (i), bad_index ());
1816             return value_data_ [it - index_data_.begin ()];
1817         }
1818
1819     public:
1820         class const_iterator;
1821         class iterator;
1822
1823         // Element lookup
1824         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1825         const_iterator find (size_type i) const {
1826             sort ();
1827             return const_iterator (*this, detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1828         }
1829         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1830         iterator find (size_type i) {
1831             sort ();
1832             return iterator (*this, detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1833         }
1834
1835
1836         class const_iterator:
1837             public container_const_reference<coordinate_vector>,
1838             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
1839                                                const_iterator, value_type> {
1840         public:
1841             typedef typename coordinate_vector::value_type value_type;
1842             typedef typename coordinate_vector::difference_type difference_type;
1843             typedef typename coordinate_vector::const_reference reference;
1844             typedef const typename coordinate_vector::pointer pointer;
1845
1846             // Construction and destruction
1847             BOOST_UBLAS_INLINE
1848             const_iterator ():
1849                 container_const_reference<self_type> (), it_ () {}
1850             BOOST_UBLAS_INLINE
1851             const_iterator (const self_type &v, const const_subiterator_type &it):
1852                 container_const_reference<self_type> (v), it_ (it) {}
1853             BOOST_UBLAS_INLINE
1854             const_iterator (const typename self_type::iterator &it):  // ISSUE self_type:: stops VC8 using std::iterator here
1855                 container_const_reference<self_type> (it ()), it_ (it.it_) {}
1856
1857             // Arithmetic
1858             BOOST_UBLAS_INLINE
1859             const_iterator &operator ++ () {
1860                 ++ it_;
1861                 return *this;
1862             }
1863             BOOST_UBLAS_INLINE
1864             const_iterator &operator -- () {
1865                 -- it_;
1866                 return *this;
1867             }
1868
1869             // Dereference
1870             BOOST_UBLAS_INLINE
1871             const_reference operator * () const {
1872                 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
1873                 return (*this) ().value_data_ [it_ - (*this) ().index_data_.begin ()];
1874             }
1875
1876             // Index
1877             BOOST_UBLAS_INLINE
1878             size_type index () const {
1879                 BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
1880                 BOOST_UBLAS_CHECK ((*this) ().zero_based (*it_) < (*this) ().size (), bad_index ());
1881                 return (*this) ().zero_based (*it_);
1882             }
1883
1884             // Assignment
1885             BOOST_UBLAS_INLINE
1886             const_iterator &operator = (const const_iterator &it) {
1887                 container_const_reference<self_type>::assign (&it ());
1888                 it_ = it.it_;
1889                 return *this;
1890             }
1891
1892             // Comparison
1893             BOOST_UBLAS_INLINE
1894             bool operator == (const const_iterator &it) const {
1895                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1896                 return it_ == it.it_;
1897             }
1898
1899         private:
1900             const_subiterator_type it_;
1901         };
1902
1903         BOOST_UBLAS_INLINE
1904         const_iterator begin () const {
1905             return find (0);
1906         }
1907         BOOST_UBLAS_INLINE
1908         const_iterator end () const {
1909             return find (size_);
1910         }
1911
1912         class iterator:
1913             public container_reference<coordinate_vector>,
1914             public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
1915                                                iterator, value_type> {
1916         public:
1917             typedef typename coordinate_vector::value_type value_type;
1918             typedef typename coordinate_vector::difference_type difference_type;
1919             typedef typename coordinate_vector::true_reference reference;
1920             typedef typename coordinate_vector::pointer pointer;
1921
1922             // Construction and destruction
1923             BOOST_UBLAS_INLINE
1924             iterator ():
1925                 container_reference<self_type> (), it_ () {}
1926             BOOST_UBLAS_INLINE
1927             iterator (self_type &v, const subiterator_type &it):
1928                 container_reference<self_type> (v), it_ (it) {}
1929
1930             // Arithmetic
1931             BOOST_UBLAS_INLINE
1932             iterator &operator ++ () {
1933                 ++ it_;
1934                 return *this;
1935             }
1936             BOOST_UBLAS_INLINE
1937             iterator &operator -- () {
1938                 -- it_;
1939                 return *this;
1940             }
1941
1942             // Dereference
1943             BOOST_UBLAS_INLINE
1944             reference operator * () const {
1945                 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
1946                 return (*this) ().value_data_ [it_ - (*this) ().index_data_.begin ()];
1947             }
1948
1949             // Index
1950             BOOST_UBLAS_INLINE
1951             size_type index () const {
1952                 BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
1953                 BOOST_UBLAS_CHECK ((*this) ().zero_based (*it_) < (*this) ().size (), bad_index ());
1954                 return (*this) ().zero_based (*it_);
1955             }
1956
1957             // Assignment
1958             BOOST_UBLAS_INLINE
1959             iterator &operator = (const iterator &it) {
1960                 container_reference<self_type>::assign (&it ());
1961                 it_ = it.it_;
1962                 return *this;
1963             }
1964
1965             // Comparison
1966             BOOST_UBLAS_INLINE
1967             bool operator == (const iterator &it) const {
1968                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1969                 return it_ == it.it_;
1970             }
1971
1972         private:
1973             subiterator_type it_;
1974
1975             friend class const_iterator;
1976         };
1977
1978         BOOST_UBLAS_INLINE
1979         iterator begin () {
1980             return find (0);
1981         }
1982         BOOST_UBLAS_INLINE
1983         iterator end () {
1984             return find (size_);
1985         }
1986
1987         // Reverse iterator
1988         typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
1989         typedef reverse_iterator_base<iterator> reverse_iterator;
1990
1991         BOOST_UBLAS_INLINE
1992         const_reverse_iterator rbegin () const {
1993             return const_reverse_iterator (end ());
1994         }
1995         BOOST_UBLAS_INLINE
1996         const_reverse_iterator rend () const {
1997             return const_reverse_iterator (begin ());
1998         }
1999         BOOST_UBLAS_INLINE
2000         reverse_iterator rbegin () {
2001             return reverse_iterator (end ());
2002         }
2003         BOOST_UBLAS_INLINE
2004         reverse_iterator rend () {
2005             return reverse_iterator (begin ());
2006         }
2007
2008          // Serialization
2009         template<class Archive>
2010         void serialize(Archive & ar, const unsigned int /* file_version */){
2011             serialization::collection_size_type s (size_);
2012             ar & serialization::make_nvp("size",s);
2013             if (Archive::is_loading::value) {
2014                 size_ = s;
2015             }
2016             // ISSUE: filled may be much less than capacity
2017             // ISSUE: index_data_ and value_data_ are undefined between filled and capacity (trouble with 'nan'-values)
2018             ar & serialization::make_nvp("capacity", capacity_);
2019             ar & serialization::make_nvp("filled", filled_);
2020             ar & serialization::make_nvp("sorted_filled", sorted_filled_);
2021             ar & serialization::make_nvp("sorted", sorted_);
2022             ar & serialization::make_nvp("index_data", index_data_);
2023             ar & serialization::make_nvp("value_data", value_data_);
2024             storage_invariants();
2025         }
2026
2027     private:
2028         void storage_invariants () const
2029         {
2030             BOOST_UBLAS_CHECK (capacity_ == index_data_.size (), internal_logic ());
2031             BOOST_UBLAS_CHECK (capacity_ == value_data_.size (), internal_logic ());
2032             BOOST_UBLAS_CHECK (filled_ <= capacity_, internal_logic ());
2033             BOOST_UBLAS_CHECK (sorted_filled_ <= filled_, internal_logic ());
2034             BOOST_UBLAS_CHECK (sorted_ == (sorted_filled_ == filled_), internal_logic ());
2035             BOOST_UBLAS_CHECK ((0 == filled_) || (zero_based(index_data_[filled_ - 1]) < size_), internal_logic ());
2036         }
2037
2038         size_type size_;
2039         size_type capacity_;
2040         mutable typename index_array_type::size_type filled_;
2041         mutable typename index_array_type::size_type sorted_filled_;
2042         mutable bool sorted_;
2043         mutable index_array_type index_data_;
2044         mutable value_array_type value_data_;
2045         static const value_type zero_;
2046
2047         BOOST_UBLAS_INLINE
2048         static size_type zero_based (size_type k_based_index) {
2049             return k_based_index - IB;
2050         }
2051         BOOST_UBLAS_INLINE
2052         static size_type k_based (size_type zero_based_index) {
2053             return zero_based_index + IB;
2054         }
2055
2056         friend class iterator;
2057         friend class const_iterator;
2058     };
2059
2060     template<class T, std::size_t IB, class IA, class TA>
2061     const typename coordinate_vector<T, IB, IA, TA>::value_type coordinate_vector<T, IB, IA, TA>::zero_ = value_type/*zero*/();
2062
2063 }}}
2064
2065 #endif