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