2 // Copyright (c) 2000-2002
3 // Joerg Walter, Mathias Koch
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)
9 // The authors gratefully acknowledge the support of
10 // GeNeSys mbH & Co. KG in producing this work.
13 #ifndef _BOOST_UBLAS_VECTOR_SPARSE_
14 #define _BOOST_UBLAS_VECTOR_SPARSE_
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>
23 // Iterators based on ideas of Jeremy Siek
25 namespace boost { namespace numeric { namespace ublas {
27 #ifdef BOOST_UBLAS_STRICT_VECTOR_SPARSE
30 class sparse_vector_element:
31 public container_reference<V> {
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;
40 // Proxied element operations
42 pointer p = (*this) ().find_element (i_);
46 d_ = value_type/*zero*/();
49 void set (const value_type &s) const {
50 pointer p = (*this) ().find_element (i_);
52 (*this) ().insert_element (i_, s);
58 // Construction and destruction
59 sparse_vector_element (vector_type &v, size_type i):
60 container_reference<vector_type> (v), i_ (i) {
63 sparse_vector_element (const sparse_vector_element &p):
64 container_reference<vector_type> (p), i_ (p.i_) {}
66 ~sparse_vector_element () {
71 sparse_vector_element &operator = (const sparse_vector_element &p) {
72 // Overide the implict copy assignment
79 sparse_vector_element &operator = (const D &d) {
85 sparse_vector_element &operator += (const D &d) {
93 sparse_vector_element &operator -= (const D &d) {
101 sparse_vector_element &operator *= (const D &d) {
109 sparse_vector_element &operator /= (const D &d) {
119 bool operator == (const D &d) const {
125 bool operator != (const D &d) const {
130 // Conversion - weak link in proxy as d_ is not a perfect alias for the element
132 operator const_reference () const {
137 // Conversion to reference - may be invalidated
139 value_type& ref () const {
140 const pointer p = (*this) ().find_element (i_);
142 return (*this) ().insert_element (i_, value_type/*zero*/());
149 mutable value_type d_;
153 * Generalise explicit reference access
157 struct element_reference {
158 typedef R& reference;
159 static reference get_reference (reference r)
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)
174 typename detail::element_reference<VER>::reference ref (VER& ver) {
175 return detail::element_reference<VER>::get_reference (ver);
178 typename detail::element_reference<VER>::reference ref (const VER& ver) {
179 return detail::element_reference<VER>::get_reference (ver);
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;
193 static const unsigned plus_complexity = type_traits<element_type>::plus_complexity;
194 static const unsigned multiplies_complexity = type_traits<element_type>::multiplies_complexity;
198 real_type real (const_reference t) {
199 return type_traits<element_type>::real (t);
203 real_type imag (const_reference t) {
204 return type_traits<element_type>::imag (t);
208 value_type conj (const_reference t) {
209 return type_traits<element_type>::conj (t);
214 real_type type_abs (const_reference t) {
215 return type_traits<element_type>::type_abs (t);
219 value_type type_sqrt (const_reference t) {
220 return type_traits<element_type>::type_sqrt (t);
225 real_type norm_1 (const_reference t) {
226 return type_traits<element_type>::norm_1 (t);
230 real_type norm_2 (const_reference t) {
231 return type_traits<element_type>::norm_2 (t);
235 real_type norm_inf (const_reference t) {
236 return type_traits<element_type>::norm_inf (t);
241 bool equals (const_reference t1, const_reference t2) {
242 return type_traits<element_type>::equals (t1, t2);
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;
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;
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;
263 /** \brief Index map based sparse vector
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.
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$.
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>.
276 * \tparam T the type of object stored in the vector (like double, float, complex, etc...)
277 * \tparam A the type of Storage array
279 template<class T, class A>
281 public vector_container<mapped_vector<T, A> > {
283 typedef T &true_reference;
285 typedef const T *const_pointer;
286 typedef mapped_vector<T, A> self_type;
288 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
289 using vector_container<self_type>::operator ();
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;
299 typedef sparse_vector_element<self_type> reference;
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;
306 // Construction and destruction
309 vector_container<self_type> (),
310 size_ (0), data_ () {}
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));
318 mapped_vector (const mapped_vector &v):
319 vector_container<self_type> (),
320 size_ (v.size_), data_ (v.data_) {}
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);
332 size_type size () const {
336 size_type nnz_capacity () const {
337 return detail::map_capacity (data ());
340 size_type nnz () const {
341 return data (). size ();
346 const array_type &data () const {
350 array_type &data () {
357 size_type restrict_capacity (size_type non_zeros) const {
358 non_zeros = (std::min) (non_zeros, size_);
363 void resize (size_type size, bool preserve = true) {
366 data ().erase (data ().lower_bound(size_), data ().end());
375 void reserve (size_type non_zeros = 0, bool preserve = true) {
376 detail::map_reserve (data (), restrict_capacity (non_zeros));
381 pointer find_element (size_type i) {
382 return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i));
385 const_pointer find_element (size_type i) const {
386 const_subiterator_type it (data ().find (i));
387 if (it == data ().end ())
389 BOOST_UBLAS_CHECK ((*it).first == i, internal_logic ()); // broken map
390 return &(*it).second;
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 ())
400 BOOST_UBLAS_CHECK ((*it).first == i, internal_logic ()); // broken map
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;
411 reference operator () (size_type i) {
412 #ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
415 BOOST_UBLAS_CHECK (i < size_, bad_index ());
416 return reference (*this, i);
421 const_reference operator [] (size_type i) const {
425 reference operator [] (size_type i) {
429 // Element assignment
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;
440 void erase_element (size_type i) {
441 subiterator_type it = data ().find (i);
442 if (it == data ().end ())
455 mapped_vector &operator = (const mapped_vector &v) {
462 template<class C> // Container assignment without temporary
464 mapped_vector &operator = (const vector_container<C> &v) {
465 resize (v ().size (), false);
470 mapped_vector &assign_temporary (mapped_vector &v) {
476 mapped_vector &operator = (const vector_expression<AE> &ae) {
477 self_type temporary (ae, detail::map_capacity (data()));
478 return assign_temporary (temporary);
482 mapped_vector &assign (const vector_expression<AE> &ae) {
483 vector_assign<scalar_assign> (*this, ae);
487 // Computed assignment
490 mapped_vector &operator += (const vector_expression<AE> &ae) {
491 self_type temporary (*this + ae, detail::map_capacity (data()));
492 return assign_temporary (temporary);
494 template<class C> // Container assignment without temporary
496 mapped_vector &operator += (const vector_container<C> &v) {
502 mapped_vector &plus_assign (const vector_expression<AE> &ae) {
503 vector_assign<scalar_plus_assign> (*this, ae);
508 mapped_vector &operator -= (const vector_expression<AE> &ae) {
509 self_type temporary (*this - ae, detail::map_capacity (data()));
510 return assign_temporary (temporary);
512 template<class C> // Container assignment without temporary
514 mapped_vector &operator -= (const vector_container<C> &v) {
520 mapped_vector &minus_assign (const vector_expression<AE> &ae) {
521 vector_assign<scalar_minus_assign> (*this, ae);
526 mapped_vector &operator *= (const AT &at) {
527 vector_assign_scalar<scalar_multiplies_assign> (*this, at);
532 mapped_vector &operator /= (const AT &at) {
533 vector_assign_scalar<scalar_divides_assign> (*this, at);
539 void swap (mapped_vector &v) {
541 std::swap (size_, v.size_);
542 data ().swap (v.data ());
546 friend void swap (mapped_vector &v1, mapped_vector &v2) {
552 // Use storage iterator
553 typedef typename A::const_iterator const_subiterator_type;
554 typedef typename A::iterator subiterator_type;
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
566 class const_iterator;
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));
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));
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> {
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;
590 // Construction and destruction
593 container_const_reference<self_type> (), it_ () {}
595 const_iterator (const self_type &v, const const_subiterator_type &it):
596 container_const_reference<self_type> (v), it_ (it) {}
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_) {}
603 const_iterator &operator ++ () {
608 const_iterator &operator -- () {
615 const_reference operator * () const {
616 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
617 return (*it_).second;
622 size_type index () const {
623 BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
624 BOOST_UBLAS_CHECK ((*it_).first < (*this) ().size (), bad_index ());
630 const_iterator &operator = (const const_iterator &it) {
631 container_const_reference<self_type>::assign (&it ());
638 bool operator == (const const_iterator &it) const {
639 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
640 return it_ == it.it_;
644 const_subiterator_type it_;
648 const_iterator begin () const {
649 return const_iterator (*this, data ().begin ());
652 const_iterator cbegin () const {
656 const_iterator end () const {
657 return const_iterator (*this, data ().end ());
660 const_iterator cend () const {
665 public container_reference<mapped_vector>,
666 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
667 iterator, value_type> {
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;
674 // Construction and destruction
677 container_reference<self_type> (), it_ () {}
679 iterator (self_type &v, const subiterator_type &it):
680 container_reference<self_type> (v), it_ (it) {}
684 iterator &operator ++ () {
689 iterator &operator -- () {
696 reference operator * () const {
697 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
698 return (*it_).second;
703 size_type index () const {
704 BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
705 BOOST_UBLAS_CHECK ((*it_).first < (*this) ().size (), bad_index ());
711 iterator &operator = (const iterator &it) {
712 container_reference<self_type>::assign (&it ());
719 bool operator == (const iterator &it) const {
720 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
721 return it_ == it.it_;
725 subiterator_type it_;
727 friend class const_iterator;
732 return iterator (*this, data ().begin ());
736 return iterator (*this, data ().end ());
740 typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
741 typedef reverse_iterator_base<iterator> reverse_iterator;
744 const_reverse_iterator rbegin () const {
745 return const_reverse_iterator (end ());
748 const_reverse_iterator crbegin () const {
752 const_reverse_iterator rend () const {
753 return const_reverse_iterator (begin ());
756 const_reverse_iterator crend () const {
760 reverse_iterator rbegin () {
761 return reverse_iterator (end ());
764 reverse_iterator rend () {
765 return reverse_iterator (begin ());
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) {
776 ar & serialization::make_nvp("data", data_);
782 static const value_type zero_;
785 template<class T, class A>
786 const typename mapped_vector<T, A>::value_type mapped_vector<T, A>::zero_ = value_type/*zero*/();
789 // Thanks to Kresimir Fresl for extending this to cover different index bases.
791 /** \brief Compressed array based sparse vector
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).
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$.
804 * Supported parameters for the adapted array (indices and values) are \c unbounded_array<> ,
805 * \c bounded_array<> and \c std::vector<>.
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>
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> > {
816 typedef T &true_reference;
818 typedef const T *const_pointer;
819 typedef compressed_vector<T, IB, IA, TA> self_type;
821 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
822 using vector_container<self_type>::operator ();
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;
833 typedef sparse_vector_element<self_type> reference;
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;
842 // Construction and destruction
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 ();
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 ();
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 ();
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);
876 size_type size () const {
880 size_type nnz_capacity () const {
884 size_type nnz () const {
890 static size_type index_base () {
894 typename index_array_type::size_type filled () const {
898 const index_array_type &index_data () const {
902 const value_array_type &value_data () const {
906 void set_filled (const typename index_array_type::size_type & filled) {
908 storage_invariants ();
911 index_array_type &index_data () {
915 value_array_type &value_data () {
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_);
929 void resize (size_type size, bool preserve = true) {
931 capacity_ = restrict_capacity (capacity_);
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)) {
941 index_data_. resize (capacity_);
942 value_data_. resize (capacity_);
945 storage_invariants ();
950 void reserve (size_type non_zeros, bool preserve = true) {
951 capacity_ = restrict_capacity (non_zeros);
953 index_data_. resize (capacity_, size_type ());
954 value_data_. resize (capacity_, value_type ());
955 filled_ = (std::min) (capacity_, filled_);
958 index_data_. resize (capacity_);
959 value_data_. resize (capacity_);
962 storage_invariants ();
967 pointer find_element (size_type i) {
968 return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i));
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))
975 return &value_data_ [it - index_data_.begin ()];
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))
985 return value_data_ [it - index_data_.begin ()];
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*/());
994 return value_data_ [it - index_data_.begin ()];
997 reference operator () (size_type i) {
998 #ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
1001 BOOST_UBLAS_CHECK (i < size_, bad_index ());
1002 return reference (*this, i);
1007 const_reference operator [] (size_type i) const {
1011 reference operator [] (size_type i) {
1015 // Element assignment
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
1026 it = index_data_.begin () + n;
1027 std::copy_backward (it, index_data_.begin () + filled_ - 1, index_data_.begin () + filled_);
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_);
1032 storage_invariants ();
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);
1045 storage_invariants ();
1052 storage_invariants ();
1057 compressed_vector &operator = (const compressed_vector &v) {
1060 capacity_ = v.capacity_;
1061 filled_ = v.filled_;
1062 index_data_ = v.index_data_;
1063 value_data_ = v.value_data_;
1065 storage_invariants ();
1068 template<class C> // Container assignment without temporary
1070 compressed_vector &operator = (const vector_container<C> &v) {
1071 resize (v ().size (), false);
1076 compressed_vector &assign_temporary (compressed_vector &v) {
1082 compressed_vector &operator = (const vector_expression<AE> &ae) {
1083 self_type temporary (ae, capacity_);
1084 return assign_temporary (temporary);
1088 compressed_vector &assign (const vector_expression<AE> &ae) {
1089 vector_assign<scalar_assign> (*this, ae);
1093 // Computed assignment
1096 compressed_vector &operator += (const vector_expression<AE> &ae) {
1097 self_type temporary (*this + ae, capacity_);
1098 return assign_temporary (temporary);
1100 template<class C> // Container assignment without temporary
1102 compressed_vector &operator += (const vector_container<C> &v) {
1108 compressed_vector &plus_assign (const vector_expression<AE> &ae) {
1109 vector_assign<scalar_plus_assign> (*this, ae);
1114 compressed_vector &operator -= (const vector_expression<AE> &ae) {
1115 self_type temporary (*this - ae, capacity_);
1116 return assign_temporary (temporary);
1118 template<class C> // Container assignment without temporary
1120 compressed_vector &operator -= (const vector_container<C> &v) {
1126 compressed_vector &minus_assign (const vector_expression<AE> &ae) {
1127 vector_assign<scalar_minus_assign> (*this, ae);
1132 compressed_vector &operator *= (const AT &at) {
1133 vector_assign_scalar<scalar_multiplies_assign> (*this, at);
1138 compressed_vector &operator /= (const AT &at) {
1139 vector_assign_scalar<scalar_divides_assign> (*this, at);
1145 void swap (compressed_vector &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_);
1153 storage_invariants ();
1156 friend void swap (compressed_vector &v1, compressed_vector &v2) {
1160 // Back element insertion and erasure
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;
1170 storage_invariants ();
1174 BOOST_UBLAS_CHECK (filled_ > 0, external_logic ());
1176 storage_invariants ();
1181 // Use index array iterator
1182 typedef typename IA::const_iterator const_subiterator_type;
1183 typedef typename IA::iterator subiterator_type;
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 ()];
1194 class const_iterator;
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> ()));
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> ()));
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> {
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;
1218 // Construction and destruction
1221 container_const_reference<self_type> (), it_ () {}
1223 const_iterator (const self_type &v, const const_subiterator_type &it):
1224 container_const_reference<self_type> (v), it_ (it) {}
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_) {}
1231 const_iterator &operator ++ () {
1236 const_iterator &operator -- () {
1243 const_reference operator * () const {
1244 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
1245 return (*this) ().value_data_ [it_ - (*this) ().index_data_.begin ()];
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_);
1258 const_iterator &operator = (const const_iterator &it) {
1259 container_const_reference<self_type>::assign (&it ());
1266 bool operator == (const const_iterator &it) const {
1267 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1268 return it_ == it.it_;
1272 const_subiterator_type it_;
1276 const_iterator begin () const {
1280 const_iterator cbegin () const {
1284 const_iterator end () const {
1285 return find (size_);
1288 const_iterator cend () const {
1293 public container_reference<compressed_vector>,
1294 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
1295 iterator, value_type> {
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;
1302 // Construction and destruction
1305 container_reference<self_type> (), it_ () {}
1307 iterator (self_type &v, const subiterator_type &it):
1308 container_reference<self_type> (v), it_ (it) {}
1312 iterator &operator ++ () {
1317 iterator &operator -- () {
1324 reference operator * () const {
1325 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
1326 return (*this) ().value_data_ [it_ - (*this) ().index_data_.begin ()];
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_);
1339 iterator &operator = (const iterator &it) {
1340 container_reference<self_type>::assign (&it ());
1347 bool operator == (const iterator &it) const {
1348 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1349 return it_ == it.it_;
1353 subiterator_type it_;
1355 friend class const_iterator;
1364 return find (size_);
1368 typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
1369 typedef reverse_iterator_base<iterator> reverse_iterator;
1372 const_reverse_iterator rbegin () const {
1373 return const_reverse_iterator (end ());
1376 const_reverse_iterator crbegin () const {
1380 const_reverse_iterator rend () const {
1381 return const_reverse_iterator (begin ());
1384 const_reverse_iterator crend () const {
1388 reverse_iterator rbegin () {
1389 return reverse_iterator (end ());
1392 reverse_iterator rend () {
1393 return reverse_iterator (begin ());
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) {
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();
1414 void storage_invariants () const
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 ());
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_;
1430 static size_type zero_based (size_type k_based_index) {
1431 return k_based_index - IB;
1434 static size_type k_based (size_type zero_based_index) {
1435 return zero_based_index + IB;
1438 friend class iterator;
1439 friend class const_iterator;
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*/();
1445 // Thanks to Kresimir Fresl for extending this to cover different index bases.
1447 /** \brief Coordimate array based sparse vector
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.
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$.
1460 * Supported parameters for the adapted array (indices and values) are \c unbounded_array<> ,
1461 * \c bounded_array<> and \c std::vector<>.
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>
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> > {
1472 typedef T &true_reference;
1474 typedef const T *const_pointer;
1475 typedef coordinate_vector<T, IB, IA, TA> self_type;
1477 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
1478 using vector_container<self_type>::operator ();
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;
1489 typedef sparse_vector_element<self_type> reference;
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;
1498 // Construction and destruction
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 ();
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 ();
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 ();
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);
1536 size_type size () const {
1540 size_type nnz_capacity () const {
1544 size_type nnz () const {
1548 // Storage accessors
1550 static size_type index_base () {
1554 typename index_array_type::size_type filled () const {
1558 const index_array_type &index_data () const {
1562 const value_array_type &value_data () const {
1566 void set_filled (const typename index_array_type::size_type &sorted, const typename index_array_type::size_type &filled) {
1567 sorted_filled_ = sorted;
1569 storage_invariants ();
1572 index_array_type &index_data () {
1576 value_array_type &value_data () {
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
1591 void resize (size_type size, bool preserve = true) {
1593 sort (); // remove duplicate elements.
1595 capacity_ = restrict_capacity (capacity_);
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)) {
1605 index_data_. resize (capacity_);
1606 value_data_. resize (capacity_);
1609 sorted_filled_ = filled_;
1610 storage_invariants ();
1614 void reserve (size_type non_zeros, bool preserve = true) {
1616 sort (); // remove duplicate elements.
1617 capacity_ = restrict_capacity (non_zeros);
1619 index_data_. resize (capacity_, size_type ());
1620 value_data_. resize (capacity_, value_type ());
1621 filled_ = (std::min) (capacity_, filled_);
1624 index_data_. resize (capacity_);
1625 value_data_. resize (capacity_);
1628 sorted_filled_ = filled_;
1629 storage_invariants ();
1634 pointer find_element (size_type i) {
1635 return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i));
1638 const_pointer find_element (size_type i) const {
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))
1643 return &value_data_ [it - index_data_.begin ()];
1648 const_reference operator () (size_type i) const {
1649 BOOST_UBLAS_CHECK (i < size_, bad_index ());
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))
1654 return value_data_ [it - index_data_.begin ()];
1657 true_reference ref (size_type i) {
1658 BOOST_UBLAS_CHECK (i < size_, bad_index ());
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*/());
1664 return value_data_ [it - index_data_.begin ()];
1667 reference operator () (size_type i) {
1668 #ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
1671 BOOST_UBLAS_CHECK (i < size_, bad_index ());
1672 return reference (*this, i);
1677 const_reference operator [] (size_type i) const {
1681 reference operator [] (size_type i) {
1685 // Element assignment
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;
1695 storage_invariants ();
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];
1704 void erase_element (size_type i) {
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);
1713 sorted_filled_ = filled_;
1715 storage_invariants ();
1722 sorted_filled_ = filled_;
1724 storage_invariants ();
1729 coordinate_vector &operator = (const coordinate_vector &v) {
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_;
1739 storage_invariants ();
1742 template<class C> // Container assignment without temporary
1744 coordinate_vector &operator = (const vector_container<C> &v) {
1745 resize (v ().size (), false);
1750 coordinate_vector &assign_temporary (coordinate_vector &v) {
1756 coordinate_vector &operator = (const vector_expression<AE> &ae) {
1757 self_type temporary (ae, capacity_);
1758 return assign_temporary (temporary);
1762 coordinate_vector &assign (const vector_expression<AE> &ae) {
1763 vector_assign<scalar_assign> (*this, ae);
1767 // Computed assignment
1770 coordinate_vector &operator += (const vector_expression<AE> &ae) {
1771 self_type temporary (*this + ae, capacity_);
1772 return assign_temporary (temporary);
1774 template<class C> // Container assignment without temporary
1776 coordinate_vector &operator += (const vector_container<C> &v) {
1782 coordinate_vector &plus_assign (const vector_expression<AE> &ae) {
1783 vector_assign<scalar_plus_assign> (*this, ae);
1788 coordinate_vector &operator -= (const vector_expression<AE> &ae) {
1789 self_type temporary (*this - ae, capacity_);
1790 return assign_temporary (temporary);
1792 template<class C> // Container assignment without temporary
1794 coordinate_vector &operator -= (const vector_container<C> &v) {
1800 coordinate_vector &minus_assign (const vector_expression<AE> &ae) {
1801 vector_assign<scalar_minus_assign> (*this, ae);
1806 coordinate_vector &operator *= (const AT &at) {
1807 vector_assign_scalar<scalar_multiplies_assign> (*this, at);
1812 coordinate_vector &operator /= (const AT &at) {
1813 vector_assign_scalar<scalar_divides_assign> (*this, at);
1819 void swap (coordinate_vector &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_);
1829 storage_invariants ();
1832 friend void swap (coordinate_vector &v1, coordinate_vector &v2) {
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 {
1839 size_type mid = (beg + end) / 2;
1840 if (index_data_[mid] < index_data_[target]) {
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;
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]);
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);
1866 rig_mid = (mid + end) / 2;
1867 lef_mid = lower_bound(beg, mid, rig_mid);
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);
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);
1878 // Sorting and summation of duplicates
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_);
1890 const typename array_pair::iterator iunsorted = ipa.begin ();
1891 std::sort (iunsorted, ipa.end ());
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]) {
1900 index_data_ [filled] = index_data_ [i];
1901 value_data_ [filled] = value_data_ [i];
1904 value_data_ [filled] += value_data_ [i];
1907 filled_ = filled + 1;
1908 sorted_filled_ = filled_;
1910 storage_invariants ();
1914 // Back element insertion and erasure
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;
1925 sorted_filled_ = filled_;
1926 storage_invariants ();
1930 // ISSUE invariants could be simpilfied if sorted required as precondition
1931 BOOST_UBLAS_CHECK (filled_ > 0, external_logic ());
1933 sorted_filled_ = (std::min) (sorted_filled_, filled_);
1934 sorted_ = sorted_filled_ = filled_;
1935 storage_invariants ();
1940 // Use index array iterator
1941 typedef typename IA::const_iterator const_subiterator_type;
1942 typedef typename IA::iterator subiterator_type;
1945 true_reference at_element (size_type i) {
1946 BOOST_UBLAS_CHECK (i < size_, bad_index ());
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 ()];
1954 class const_iterator;
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 {
1961 return const_iterator (*this, detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
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) {
1966 return iterator (*this, detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
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> {
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;
1980 // Construction and destruction
1983 container_const_reference<self_type> (), it_ () {}
1985 const_iterator (const self_type &v, const const_subiterator_type &it):
1986 container_const_reference<self_type> (v), it_ (it) {}
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_) {}
1993 const_iterator &operator ++ () {
1998 const_iterator &operator -- () {
2005 const_reference operator * () const {
2006 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
2007 return (*this) ().value_data_ [it_ - (*this) ().index_data_.begin ()];
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_);
2020 const_iterator &operator = (const const_iterator &it) {
2021 container_const_reference<self_type>::assign (&it ());
2028 bool operator == (const const_iterator &it) const {
2029 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2030 return it_ == it.it_;
2034 const_subiterator_type it_;
2038 const_iterator begin () const {
2042 const_iterator cbegin () const {
2046 const_iterator end () const {
2047 return find (size_);
2050 const_iterator cend () const {
2055 public container_reference<coordinate_vector>,
2056 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
2057 iterator, value_type> {
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;
2064 // Construction and destruction
2067 container_reference<self_type> (), it_ () {}
2069 iterator (self_type &v, const subiterator_type &it):
2070 container_reference<self_type> (v), it_ (it) {}
2074 iterator &operator ++ () {
2079 iterator &operator -- () {
2086 reference operator * () const {
2087 BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
2088 return (*this) ().value_data_ [it_ - (*this) ().index_data_.begin ()];
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_);
2101 iterator &operator = (const iterator &it) {
2102 container_reference<self_type>::assign (&it ());
2109 bool operator == (const iterator &it) const {
2110 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2111 return it_ == it.it_;
2115 subiterator_type it_;
2117 friend class const_iterator;
2126 return find (size_);
2130 typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
2131 typedef reverse_iterator_base<iterator> reverse_iterator;
2134 const_reverse_iterator rbegin () const {
2135 return const_reverse_iterator (end ());
2138 const_reverse_iterator crbegin () const {
2142 const_reverse_iterator rend () const {
2143 return const_reverse_iterator (begin ());
2146 const_reverse_iterator crend () const {
2150 reverse_iterator rbegin () {
2151 return reverse_iterator (end ());
2154 reverse_iterator rend () {
2155 return reverse_iterator (begin ());
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) {
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();
2178 void storage_invariants () const
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 ());
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_;
2198 static size_type zero_based (size_type k_based_index) {
2199 return k_based_index - IB;
2202 static size_type k_based (size_type zero_based_index) {
2203 return zero_based_index + IB;
2206 friend class iterator;
2207 friend class const_iterator;
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*/();