1	//
2	// Copyright (c) 2000-2007
3	// Joerg Walter, Mathias Koch, Gunter Winkler
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_MATRIX_SPARSE_
14	#define _BOOST_UBLAS_MATRIX_SPARSE_
15
16	#include <boost/numeric/ublas/vector_sparse.hpp>
17	#include <boost/numeric/ublas/matrix_expression.hpp>
18	#include <boost/numeric/ublas/detail/matrix_assign.hpp>
19	#if BOOST_UBLAS_TYPE_CHECK
20	#include <boost/numeric/ublas/matrix.hpp>
21	#endif
22
23	// Iterators based on ideas of Jeremy Siek
24
25	namespace boost { namespace numeric { namespace ublas {
26
27	#ifdef BOOST_UBLAS_STRICT_MATRIX_SPARSE
28
29	template<class M>
30	class sparse_matrix_element:
31	public container_reference<M> {
32	public:
33	typedef M matrix_type;
34	typedef typename M::size_type size_type;
35	typedef typename M::value_type value_type;
36	typedef const value_type &const_reference;
37	typedef value_type *pointer;
38	typedef const value_type *const_pointer;
39
40	private:
41	// Proxied element operations
42	void get_d () const {
43	const_pointer p = (*this) ().find_element (i_, j_);
44	if (p)
45	d_ = *p;
46	else
47	d_ = value_type/*zero*/();
48	}
49
50	void set (const value_type &s) const {
51	pointer p = (*this) ().find_element (i_, j_);
52	if (!p)
53	(*this) ().insert_element (i_, j_, s);
54	else
55	*p = s;
56	}
57
58	public:
59	// Construction and destruction
60	BOOST_UBLAS_INLINE
61	sparse_matrix_element (matrix_type &m, size_type i, size_type j):
62	container_reference<matrix_type> (m), i_ (i), j_ (j) {
63	}
64	BOOST_UBLAS_INLINE
65	sparse_matrix_element (const sparse_matrix_element &p):
66	container_reference<matrix_type> (p), i_ (p.i_), j_ (p.j_) {}
67	BOOST_UBLAS_INLINE
68	~sparse_matrix_element () {
69	}
70
71	// Assignment
72	BOOST_UBLAS_INLINE
73	sparse_matrix_element &operator = (const sparse_matrix_element &p) {
74	// Overide the implict copy assignment
75	p.get_d ();
76	set (p.d_);
77	return *this;
78	}
79	template<class D>
80	BOOST_UBLAS_INLINE
81	sparse_matrix_element &operator = (const D &d) {
82	set (d);
83	return *this;
84	}
85	template<class D>
86	BOOST_UBLAS_INLINE
87	sparse_matrix_element &operator += (const D &d) {
88	get_d ();
89	d_ += d;
90	set (d_);
91	return *this;
92	}
93	template<class D>
94	BOOST_UBLAS_INLINE
95	sparse_matrix_element &operator -= (const D &d) {
96	get_d ();
97	d_ -= d;
98	set (d_);
99	return *this;
100	}
101	template<class D>
102	BOOST_UBLAS_INLINE
103	sparse_matrix_element &operator *= (const D &d) {
104	get_d ();
105	d_ *= d;
106	set (d_);
107	return *this;
108	}
109	template<class D>
110	BOOST_UBLAS_INLINE
111	sparse_matrix_element &operator /= (const D &d) {
112	get_d ();
113	d_ /= d;
114	set (d_);
115	return *this;
116	}
117
118	// Comparison
119	template<class D>
120	BOOST_UBLAS_INLINE
121	bool operator == (const D &d) const {
122	get_d ();
123	return d_ == d;
124	}
125	template<class D>
126	BOOST_UBLAS_INLINE
127	bool operator != (const D &d) const {
128	get_d ();
129	return d_ != d;
130	}
131
132	// Conversion - weak link in proxy as d_ is not a perfect alias for the element
133	BOOST_UBLAS_INLINE
134	operator const_reference () const {
135	get_d ();
136	return d_;
137	}
138
139	// Conversion to reference - may be invalidated
140	BOOST_UBLAS_INLINE
141	value_type& ref () const {
142	const pointer p = (*this) ().find_element (i_, j_);
143	if (!p)
144	return (*this) ().insert_element (i_, j_, value_type/*zero*/());
145	else
146	return *p;
147	}
148
149	private:
150	size_type i_;
151	size_type j_;
152	mutable value_type d_;
153	};
154
155	/*
156	* Generalise explicit reference access
157	*/
158	namespace detail {
159	template <class V>
160	struct element_reference<sparse_matrix_element<V> > {
161	typedef typename V::value_type& reference;
162	static reference get_reference (const sparse_matrix_element<V>& sve)
163	{
164	return sve.ref ();
165	}
166	};
167	}
168
169
170	template<class M>
171	struct type_traits<sparse_matrix_element<M> > {
172	typedef typename M::value_type element_type;
173	typedef type_traits<sparse_matrix_element<M> > self_type;
174	typedef typename type_traits<element_type>::value_type value_type;
175	typedef typename type_traits<element_type>::const_reference const_reference;
176	typedef sparse_matrix_element<M> reference;
177	typedef typename type_traits<element_type>::real_type real_type;
178	typedef typename type_traits<element_type>::precision_type precision_type;
179
180	static const unsigned plus_complexity = type_traits<element_type>::plus_complexity;
181	static const unsigned multiplies_complexity = type_traits<element_type>::multiplies_complexity;
182
183	static
184	BOOST_UBLAS_INLINE
185	real_type real (const_reference t) {
186	return type_traits<element_type>::real (t);
187	}
188	static
189	BOOST_UBLAS_INLINE
190	real_type imag (const_reference t) {
191	return type_traits<element_type>::imag (t);
192	}
193	static
194	BOOST_UBLAS_INLINE
195	value_type conj (const_reference t) {
196	return type_traits<element_type>::conj (t);
197	}
198
199	static
200	BOOST_UBLAS_INLINE
201	real_type type_abs (const_reference t) {
202	return type_traits<element_type>::type_abs (t);
203	}
204	static
205	BOOST_UBLAS_INLINE
206	value_type type_sqrt (const_reference t) {
207	return type_traits<element_type>::type_sqrt (t);
208	}
209
210	static
211	BOOST_UBLAS_INLINE
212	real_type norm_1 (const_reference t) {
213	return type_traits<element_type>::norm_1 (t);
214	}
215	static
216	BOOST_UBLAS_INLINE
217	real_type norm_2 (const_reference t) {
218	return type_traits<element_type>::norm_2 (t);
219	}
220	static
221	BOOST_UBLAS_INLINE
222	real_type norm_inf (const_reference t) {
223	return type_traits<element_type>::norm_inf (t);
224	}
225
226	static
227	BOOST_UBLAS_INLINE
228	bool equals (const_reference t1, const_reference t2) {
229	return type_traits<element_type>::equals (t1, t2);
230	}
231	};
232
233	template<class M1, class T2>
234	struct promote_traits<sparse_matrix_element<M1>, T2> {
235	typedef typename promote_traits<typename sparse_matrix_element<M1>::value_type, T2>::promote_type promote_type;
236	};
237	template<class T1, class M2>
238	struct promote_traits<T1, sparse_matrix_element<M2> > {
239	typedef typename promote_traits<T1, typename sparse_matrix_element<M2>::value_type>::promote_type promote_type;
240	};
241	template<class M1, class M2>
242	struct promote_traits<sparse_matrix_element<M1>, sparse_matrix_element<M2> > {
243	typedef typename promote_traits<typename sparse_matrix_element<M1>::value_type,
244	typename sparse_matrix_element<M2>::value_type>::promote_type promote_type;
245	};
246
247	#endif
248
249	/** \brief Index map based sparse matrix of values of type \c T
250	*
251	* This class represents a matrix by using a \c key to value mapping. The default type is
252	* \code template<class T, class L = row_major, class A = map_std<std::size_t, T> > class mapped_matrix; \endcode
253	* So, by default a STL map container is used to associate keys and values. The key is computed depending on
254	* the layout type \c L as \code key = layout_type::element(i, size1_, j, size2_); \endcode
255	* which means \code key = (i*size2+j) \endcode for a row major matrix.
256	* Limitations: The matrix size must not exceed \f$(size1*size2) < \f$ \code std::limits<std::size_t> \endcode.
257	* The \ref find1() and \ref find2() operations have a complexity of at least \f$\mathcal{O}(log(nnz))\f$, depending
258	* on the efficiency of \c std::lower_bound on the key set of the map.
259	* Orientation and storage can also be specified, otherwise a row major orientation is used.
260	* It is \b not required by the storage to initialize elements of the matrix. By default, the orientation is \c row_major.
261	*
262	* \sa fwd.hpp, storage_sparse.hpp
263	*
264	* \tparam T the type of object stored in the matrix (like double, float, complex, etc...)
265	* \tparam L the storage organization. It can be either \c row_major or \c column_major. By default it is \c row_major
266	*/
267	template<class T, class L, class A>
268	class mapped_matrix:
269	public matrix_container<mapped_matrix<T, L, A> > {
270
271	typedef T &true_reference;
272	typedef T *pointer;
273	typedef const T * const_pointer;
274	typedef L layout_type;
275	typedef mapped_matrix<T, L, A> self_type;
276	public:
277	#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
278	using matrix_container<self_type>::operator ();
279	#endif
280	typedef typename A::size_type size_type;
281	typedef typename A::difference_type difference_type;
282	typedef T value_type;
283	typedef A array_type;
284	typedef const T &const_reference;
285	#ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
286	typedef typename detail::map_traits<A, T>::reference reference;
287	#else
288	typedef sparse_matrix_element<self_type> reference;
289	#endif
290	typedef const matrix_reference<const self_type> const_closure_type;
291	typedef matrix_reference<self_type> closure_type;
292	typedef mapped_vector<T, A> vector_temporary_type;
293	typedef self_type matrix_temporary_type;
294	typedef sparse_tag storage_category;
295	typedef typename L::orientation_category orientation_category;
296
297	// Construction and destruction
298	BOOST_UBLAS_INLINE
299	mapped_matrix ():
300	matrix_container<self_type> (),
301	size1_ (0), size2_ (0), data_ () {}
302	BOOST_UBLAS_INLINE
303	mapped_matrix (size_type size1, size_type size2, size_type non_zeros = 0):
304	matrix_container<self_type> (),
305	size1_ (size1), size2_ (size2), data_ () {
306	detail::map_reserve (data (), restrict_capacity (non_zeros));
307	}
308	BOOST_UBLAS_INLINE
309	mapped_matrix (const mapped_matrix &m):
310	matrix_container<self_type> (),
311	size1_ (m.size1_), size2_ (m.size2_), data_ (m.data_) {}
312	template<class AE>
313	BOOST_UBLAS_INLINE
314	mapped_matrix (const matrix_expression<AE> &ae, size_type non_zeros = 0):
315	matrix_container<self_type> (),
316	size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), data_ () {
317	detail::map_reserve (data (), restrict_capacity (non_zeros));
318	matrix_assign<scalar_assign> (*this, ae);
319	}
320
321	// Accessors
322	BOOST_UBLAS_INLINE
323	size_type size1 () const {
324	return size1_;
325	}
326	BOOST_UBLAS_INLINE
327	size_type size2 () const {
328	return size2_;
329	}
330	BOOST_UBLAS_INLINE
331	size_type nnz_capacity () const {
332	return detail::map_capacity (data ());
333	}
334	BOOST_UBLAS_INLINE
335	size_type nnz () const {
336	return data (). size ();
337	}
338
339	// Storage accessors
340	BOOST_UBLAS_INLINE
341	const array_type &data () const {
342	return data_;
343	}
344	BOOST_UBLAS_INLINE
345	array_type &data () {
346	return data_;
347	}
348
349	// Resizing
350	private:
351	BOOST_UBLAS_INLINE
352	size_type restrict_capacity (size_type non_zeros) const {
353	// Guarding against overflow - thanks to Alexei Novakov for the hint.
354	// non_zeros = (std::min) (non_zeros, size1_ * size2_);
355	if (size1_ > 0 && non_zeros / size1_ >= size2_)
356	non_zeros = size1_ * size2_;
357	return non_zeros;
358	}
359	public:
360	BOOST_UBLAS_INLINE
361	void resize (size_type size1, size_type size2, bool preserve = true) {
362	// FIXME preserve unimplemented
363	BOOST_UBLAS_CHECK (!preserve, internal_logic ());
364	size1_ = size1;
365	size2_ = size2;
366	data ().clear ();
367	}
368
369	// Reserving
370	BOOST_UBLAS_INLINE
371	void reserve (size_type non_zeros, bool /*preserve*/ = true) {
372	detail::map_reserve (data (), restrict_capacity (non_zeros));
373	}
374
375	// Element support
376	BOOST_UBLAS_INLINE
377	pointer find_element (size_type i, size_type j) {
378	return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i, j));
379	}
380	BOOST_UBLAS_INLINE
381	const_pointer find_element (size_type i, size_type j) const {
382	const size_type element = layout_type::element (i, size1_, j, size2_);
383	const_subiterator_type it (data ().find (element));
384	if (it == data ().end ())
385	return 0;
386	BOOST_UBLAS_CHECK ((*it).first == element, internal_logic ()); // broken map
387	return &(*it).second;
388	}
389
390	// Element access
391	BOOST_UBLAS_INLINE
392	const_reference operator () (size_type i, size_type j) const {
393	const size_type element = layout_type::element (i, size1_, j, size2_);
394	const_subiterator_type it (data ().find (element));
395	if (it == data ().end ())
396	return zero_;
397	BOOST_UBLAS_CHECK ((*it).first == element, internal_logic ()); // broken map
398	return (*it).second;
399	}
400	BOOST_UBLAS_INLINE
401	reference operator () (size_type i, size_type j) {
402	#ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
403	const size_type element = layout_type::element (i, size1_, j, size2_);
404	std::pair<subiterator_type, bool> ii (data ().insert (typename array_type::value_type (element, value_type/*zero*/())));
405	BOOST_UBLAS_CHECK ((ii.first)->first == element, internal_logic ()); // broken map
406	return (ii.first)->second;
407	#else
408	return reference (*this, i, j);
409	#endif
410	}
411
412	// Element assingment
413	BOOST_UBLAS_INLINE
414	true_reference insert_element (size_type i, size_type j, const_reference t) {
415	BOOST_UBLAS_CHECK (!find_element (i, j), bad_index ()); // duplicate element
416	const size_type element = layout_type::element (i, size1_, j, size2_);
417	std::pair<subiterator_type, bool> ii (data ().insert (typename array_type::value_type (element, t)));
418	BOOST_UBLAS_CHECK ((ii.first)->first == element, internal_logic ()); // broken map
419	if (!ii.second) // existing element
420	(ii.first)->second = t;
421	return (ii.first)->second;
422	}
423	BOOST_UBLAS_INLINE
424	void erase_element (size_type i, size_type j) {
425	subiterator_type it = data ().find (layout_type::element (i, size1_, j, size2_));
426	if (it == data ().end ())
427	return;
428	data ().erase (it);
429	}
430
431	// Zeroing
432	BOOST_UBLAS_INLINE
433	void clear () {
434	data ().clear ();
435	}
436
437	// Assignment
438	BOOST_UBLAS_INLINE
439	mapped_matrix &operator = (const mapped_matrix &m) {
440	if (this != &m) {
441	size1_ = m.size1_;
442	size2_ = m.size2_;
443	data () = m.data ();
444	}
445	return *this;
446	}
447	template<class C> // Container assignment without temporary
448	BOOST_UBLAS_INLINE
449	mapped_matrix &operator = (const matrix_container<C> &m) {
450	resize (size1: m ().size1 (), size2: m ().size2 (), preserve: false);
451	assign (m);
452	return *this;
453	}
454	BOOST_UBLAS_INLINE
455	mapped_matrix &assign_temporary (mapped_matrix &m) {
456	swap (m);
457	return *this;
458	}
459	template<class AE>
460	BOOST_UBLAS_INLINE
461	mapped_matrix &operator = (const matrix_expression<AE> &ae) {
462	self_type temporary (ae, detail::map_capacity (data ()));
463	return assign_temporary (m&: temporary);
464	}
465	template<class AE>
466	BOOST_UBLAS_INLINE
467	mapped_matrix &assign (const matrix_expression<AE> &ae) {
468	matrix_assign<scalar_assign> (*this, ae);
469	return *this;
470	}
471	template<class AE>
472	BOOST_UBLAS_INLINE
473	mapped_matrix& operator += (const matrix_expression<AE> &ae) {
474	self_type temporary (*this + ae, detail::map_capacity (data ()));
475	return assign_temporary (m&: temporary);
476	}
477	template<class C> // Container assignment without temporary
478	BOOST_UBLAS_INLINE
479	mapped_matrix &operator += (const matrix_container<C> &m) {
480	plus_assign (m);
481	return *this;
482	}
483	template<class AE>
484	BOOST_UBLAS_INLINE
485	mapped_matrix &plus_assign (const matrix_expression<AE> &ae) {
486	matrix_assign<scalar_plus_assign> (*this, ae);
487	return *this;
488	}
489	template<class AE>
490	BOOST_UBLAS_INLINE
491	mapped_matrix& operator -= (const matrix_expression<AE> &ae) {
492	self_type temporary (*this - ae, detail::map_capacity (data ()));
493	return assign_temporary (m&: temporary);
494	}
495	template<class C> // Container assignment without temporary
496	BOOST_UBLAS_INLINE
497	mapped_matrix &operator -= (const matrix_container<C> &m) {
498	minus_assign (m);
499	return *this;
500	}
501	template<class AE>
502	BOOST_UBLAS_INLINE
503	mapped_matrix &minus_assign (const matrix_expression<AE> &ae) {
504	matrix_assign<scalar_minus_assign> (*this, ae);
505	return *this;
506	}
507	template<class AT>
508	BOOST_UBLAS_INLINE
509	mapped_matrix& operator *= (const AT &at) {
510	matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
511	return *this;
512	}
513	template<class AT>
514	BOOST_UBLAS_INLINE
515	mapped_matrix& operator /= (const AT &at) {
516	matrix_assign_scalar<scalar_divides_assign> (*this, at);
517	return *this;
518	}
519
520	// Swapping
521	BOOST_UBLAS_INLINE
522	void swap (mapped_matrix &m) {
523	if (this != &m) {
524	std::swap (size1_, m.size1_);
525	std::swap (size2_, m.size2_);
526	data ().swap (m.data ());
527	}
528	}
529	BOOST_UBLAS_INLINE
530	friend void swap (mapped_matrix &m1, mapped_matrix &m2) {
531	m1.swap (m2);
532	}
533
534	// Iterator types
535	private:
536	// Use storage iterator
537	typedef typename A::const_iterator const_subiterator_type;
538	typedef typename A::iterator subiterator_type;
539
540	BOOST_UBLAS_INLINE
541	true_reference at_element (size_type i, size_type j) {
542	const size_type element = layout_type::element (i, size1_, j, size2_);
543	subiterator_type it (data ().find (element));
544	BOOST_UBLAS_CHECK (it != data ().end(), bad_index ());
545	BOOST_UBLAS_CHECK ((*it).first == element, internal_logic ()); // broken map
546	return it->second;
547	}
548
549	public:
550	class const_iterator1;
551	class iterator1;
552	class const_iterator2;
553	class iterator2;
554	typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
555	typedef reverse_iterator_base1<iterator1> reverse_iterator1;
556	typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
557	typedef reverse_iterator_base2<iterator2> reverse_iterator2;
558
559	// Element lookup
560	// BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
561	const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const {
562	const_subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_)));
563	const_subiterator_type it_end (data ().end ());
564	size_type index1 = size_type (-1);
565	size_type index2 = size_type (-1);
566	while (rank == 1 && it != it_end) {
567	index1 = layout_type::index_i ((*it).first, size1_, size2_);
568	index2 = layout_type::index_j ((*it).first, size1_, size2_);
569	if (direction > 0) {
570	if ((index1 >= i && index2 == j) || (i >= size1_))
571	break;
572	++ i;
573	} else /* if (direction < 0) */ {
574	if ((index1 <= i && index2 == j) || (i == 0))
575	break;
576	-- i;
577	}
578	it = data ().lower_bound (layout_type::address (i, size1_, j, size2_));
579	}
580	if (rank == 1 && index2 != j) {
581	if (direction > 0)
582	i = size1_;
583	else /* if (direction < 0) */
584	i = 0;
585	rank = 0;
586	}
587	return const_iterator1 (*this, rank, i, j, it);
588	}
589	// BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
590	iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) {
591	subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_)));
592	subiterator_type it_end (data ().end ());
593	size_type index1 = size_type (-1);
594	size_type index2 = size_type (-1);
595	while (rank == 1 && it != it_end) {
596	index1 = layout_type::index_i ((*it).first, size1_, size2_);
597	index2 = layout_type::index_j ((*it).first, size1_, size2_);
598	if (direction > 0) {
599	if ((index1 >= i && index2 == j) || (i >= size1_))
600	break;
601	++ i;
602	} else /* if (direction < 0) */ {
603	if ((index1 <= i && index2 == j) || (i == 0))
604	break;
605	-- i;
606	}
607	it = data ().lower_bound (layout_type::address (i, size1_, j, size2_));
608	}
609	if (rank == 1 && index2 != j) {
610	if (direction > 0)
611	i = size1_;
612	else /* if (direction < 0) */
613	i = 0;
614	rank = 0;
615	}
616	return iterator1 (*this, rank, i, j, it);
617	}
618	// BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
619	const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const {
620	const_subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_)));
621	const_subiterator_type it_end (data ().end ());
622	size_type index1 = size_type (-1);
623	size_type index2 = size_type (-1);
624	while (rank == 1 && it != it_end) {
625	index1 = layout_type::index_i ((*it).first, size1_, size2_);
626	index2 = layout_type::index_j ((*it).first, size1_, size2_);
627	if (direction > 0) {
628	if ((index2 >= j && index1 == i) || (j >= size2_))
629	break;
630	++ j;
631	} else /* if (direction < 0) */ {
632	if ((index2 <= j && index1 == i) || (j == 0))
633	break;
634	-- j;
635	}
636	it = data ().lower_bound (layout_type::address (i, size1_, j, size2_));
637	}
638	if (rank == 1 && index1 != i) {
639	if (direction > 0)
640	j = size2_;
641	else /* if (direction < 0) */
642	j = 0;
643	rank = 0;
644	}
645	return const_iterator2 (*this, rank, i, j, it);
646	}
647	// BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
648	iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) {
649	subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_)));
650	subiterator_type it_end (data ().end ());
651	size_type index1 = size_type (-1);
652	size_type index2 = size_type (-1);
653	while (rank == 1 && it != it_end) {
654	index1 = layout_type::index_i ((*it).first, size1_, size2_);
655	index2 = layout_type::index_j ((*it).first, size1_, size2_);
656	if (direction > 0) {
657	if ((index2 >= j && index1 == i) || (j >= size2_))
658	break;
659	++ j;
660	} else /* if (direction < 0) */ {
661	if ((index2 <= j && index1 == i) || (j == 0))
662	break;
663	-- j;
664	}
665	it = data ().lower_bound (layout_type::address (i, size1_, j, size2_));
666	}
667	if (rank == 1 && index1 != i) {
668	if (direction > 0)
669	j = size2_;
670	else /* if (direction < 0) */
671	j = 0;
672	rank = 0;
673	}
674	return iterator2 (*this, rank, i, j, it);
675	}
676
677
678	class const_iterator1:
679	public container_const_reference<mapped_matrix>,
680	public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
681	const_iterator1, value_type> {
682	public:
683	typedef typename mapped_matrix::value_type value_type;
684	typedef typename mapped_matrix::difference_type difference_type;
685	typedef typename mapped_matrix::const_reference reference;
686	typedef const typename mapped_matrix::pointer pointer;
687
688	typedef const_iterator2 dual_iterator_type;
689	typedef const_reverse_iterator2 dual_reverse_iterator_type;
690
691	// Construction and destruction
692	BOOST_UBLAS_INLINE
693	const_iterator1 ():
694	container_const_reference<self_type> (), rank_ (), i_ (), j_ (), it_ () {}
695	BOOST_UBLAS_INLINE
696	const_iterator1 (const self_type &m, int rank, size_type i, size_type j, const const_subiterator_type &it):
697	container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), it_ (it) {}
698	BOOST_UBLAS_INLINE
699	const_iterator1 (const iterator1 &it):
700	container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), it_ (it.it_) {}
701
702	// Arithmetic
703	BOOST_UBLAS_INLINE
704	const_iterator1 &operator ++ () {
705	if (rank_ == 1 && layout_type::fast_i ())
706	++ it_;
707	else
708	*this = (*this) ().find1 (rank_, index1 () + 1, j_, 1);
709	return *this;
710	}
711	BOOST_UBLAS_INLINE
712	const_iterator1 &operator -- () {
713	if (rank_ == 1 && layout_type::fast_i ())
714	-- it_;
715	else
716	*this = (*this) ().find1 (rank_, index1 () - 1, j_, -1);
717	return *this;
718	}
719
720	// Dereference
721	BOOST_UBLAS_INLINE
722	const_reference operator * () const {
723	BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
724	BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
725	if (rank_ == 1) {
726	return (*it_).second;
727	} else {
728	return (*this) () (i_, j_);
729	}
730	}
731
732	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
733	BOOST_UBLAS_INLINE
734	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
735	typename self_type::
736	#endif
737	const_iterator2 begin () const {
738	const self_type &m = (*this) ();
739	return m.find2 (1, index1 (), 0);
740	}
741	BOOST_UBLAS_INLINE
742	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
743	typename self_type::
744	#endif
745	const_iterator2 cbegin () const {
746	return begin ();
747	}
748	BOOST_UBLAS_INLINE
749	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
750	typename self_type::
751	#endif
752	const_iterator2 end () const {
753	const self_type &m = (*this) ();
754	return m.find2 (1, index1 (), m.size2 ());
755	}
756	BOOST_UBLAS_INLINE
757	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
758	typename self_type::
759	#endif
760	const_iterator2 cend () const {
761	return end ();
762	}
763	BOOST_UBLAS_INLINE
764	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
765	typename self_type::
766	#endif
767	const_reverse_iterator2 rbegin () const {
768	return const_reverse_iterator2 (end ());
769	}
770	BOOST_UBLAS_INLINE
771	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
772	typename self_type::
773	#endif
774	const_reverse_iterator2 crbegin () const {
775	return rbegin ();
776	}
777	BOOST_UBLAS_INLINE
778	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
779	typename self_type::
780	#endif
781	const_reverse_iterator2 rend () const {
782	return const_reverse_iterator2 (begin ());
783	}
784	BOOST_UBLAS_INLINE
785	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
786	typename self_type::
787	#endif
788	const_reverse_iterator2 crend () const {
789	return rend ();
790	}
791	#endif
792
793	// Indices
794	BOOST_UBLAS_INLINE
795	size_type index1 () const {
796	BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
797	if (rank_ == 1) {
798	const self_type &m = (*this) ();
799	BOOST_UBLAS_CHECK (layout_type::index_i ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size1 (), bad_index ());
800	return layout_type::index_i ((*it_).first, m.size1 (), m.size2 ());
801	} else {
802	return i_;
803	}
804	}
805	BOOST_UBLAS_INLINE
806	size_type index2 () const {
807	if (rank_ == 1) {
808	const self_type &m = (*this) ();
809	BOOST_UBLAS_CHECK (layout_type::index_j ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size2 (), bad_index ());
810	return layout_type::index_j ((*it_).first, m.size1 (), m.size2 ());
811	} else {
812	return j_;
813	}
814	}
815
816	// Assignment
817	BOOST_UBLAS_INLINE
818	const_iterator1 &operator = (const const_iterator1 &it) {
819	container_const_reference<self_type>::assign (&it ());
820	rank_ = it.rank_;
821	i_ = it.i_;
822	j_ = it.j_;
823	it_ = it.it_;
824	return *this;
825	}
826
827	// Comparison
828	BOOST_UBLAS_INLINE
829	bool operator == (const const_iterator1 &it) const {
830	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
831	// BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
832	if (rank_ == 1 || it.rank_ == 1) {
833	return it_ == it.it_;
834	} else {
835	return i_ == it.i_ && j_ == it.j_;
836	}
837	}
838
839	private:
840	int rank_;
841	size_type i_;
842	size_type j_;
843	const_subiterator_type it_;
844	};
845
846	BOOST_UBLAS_INLINE
847	const_iterator1 begin1 () const {
848	return find1 (0, 0, 0);
849	}
850	BOOST_UBLAS_INLINE
851	const_iterator1 cbegin1 () const {
852	return begin1 ();
853	}
854	BOOST_UBLAS_INLINE
855	const_iterator1 end1 () const {
856	return find1 (0, size1_, 0);
857	}
858	BOOST_UBLAS_INLINE
859	const_iterator1 cend1 () const {
860	return end1 ();
861	}
862
863	class iterator1:
864	public container_reference<mapped_matrix>,
865	public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
866	iterator1, value_type> {
867	public:
868	typedef typename mapped_matrix::value_type value_type;
869	typedef typename mapped_matrix::difference_type difference_type;
870	typedef typename mapped_matrix::true_reference reference;
871	typedef typename mapped_matrix::pointer pointer;
872
873	typedef iterator2 dual_iterator_type;
874	typedef reverse_iterator2 dual_reverse_iterator_type;
875
876	// Construction and destruction
877	BOOST_UBLAS_INLINE
878	iterator1 ():
879	container_reference<self_type> (), rank_ (), i_ (), j_ (), it_ () {}
880	BOOST_UBLAS_INLINE
881	iterator1 (self_type &m, int rank, size_type i, size_type j, const subiterator_type &it):
882	container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), it_ (it) {}
883
884	// Arithmetic
885	BOOST_UBLAS_INLINE
886	iterator1 &operator ++ () {
887	if (rank_ == 1 && layout_type::fast_i ())
888	++ it_;
889	else
890	*this = (*this) ().find1 (rank_, index1 () + 1, j_, 1);
891	return *this;
892	}
893	BOOST_UBLAS_INLINE
894	iterator1 &operator -- () {
895	if (rank_ == 1 && layout_type::fast_i ())
896	-- it_;
897	else
898	*this = (*this) ().find1 (rank_, index1 () - 1, j_, -1);
899	return *this;
900	}
901
902	// Dereference
903	BOOST_UBLAS_INLINE
904	reference operator * () const {
905	BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
906	BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
907	if (rank_ == 1) {
908	return (*it_).second;
909	} else {
910	return (*this) ().at_element (i_, j_);
911	}
912	}
913
914	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
915	BOOST_UBLAS_INLINE
916	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
917	typename self_type::
918	#endif
919	iterator2 begin () const {
920	self_type &m = (*this) ();
921	return m.find2 (1, index1 (), 0);
922	}
923	BOOST_UBLAS_INLINE
924	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
925	typename self_type::
926	#endif
927	iterator2 end () const {
928	self_type &m = (*this) ();
929	return m.find2 (1, index1 (), m.size2 ());
930	}
931	BOOST_UBLAS_INLINE
932	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
933	typename self_type::
934	#endif
935	reverse_iterator2 rbegin () const {
936	return reverse_iterator2 (end ());
937	}
938	BOOST_UBLAS_INLINE
939	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
940	typename self_type::
941	#endif
942	reverse_iterator2 rend () const {
943	return reverse_iterator2 (begin ());
944	}
945	#endif
946
947	// Indices
948	BOOST_UBLAS_INLINE
949	size_type index1 () const {
950	BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
951	if (rank_ == 1) {
952	const self_type &m = (*this) ();
953	BOOST_UBLAS_CHECK (layout_type::index_i ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size1 (), bad_index ());
954	return layout_type::index_i ((*it_).first, m.size1 (), m.size2 ());
955	} else {
956	return i_;
957	}
958	}
959	BOOST_UBLAS_INLINE
960	size_type index2 () const {
961	if (rank_ == 1) {
962	const self_type &m = (*this) ();
963	BOOST_UBLAS_CHECK (layout_type::index_j ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size2 (), bad_index ());
964	return layout_type::index_j ((*it_).first, m.size1 (), m.size2 ());
965	} else {
966	return j_;
967	}
968	}
969
970	// Assignment
971	BOOST_UBLAS_INLINE
972	iterator1 &operator = (const iterator1 &it) {
973	container_reference<self_type>::assign (&it ());
974	rank_ = it.rank_;
975	i_ = it.i_;
976	j_ = it.j_;
977	it_ = it.it_;
978	return *this;
979	}
980
981	// Comparison
982	BOOST_UBLAS_INLINE
983	bool operator == (const iterator1 &it) const {
984	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
985	// BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
986	if (rank_ == 1 || it.rank_ == 1) {
987	return it_ == it.it_;
988	} else {
989	return i_ == it.i_ && j_ == it.j_;
990	}
991	}
992
993	private:
994	int rank_;
995	size_type i_;
996	size_type j_;
997	subiterator_type it_;
998
999	friend class const_iterator1;
1000	};
1001
1002	BOOST_UBLAS_INLINE
1003	iterator1 begin1 () {
1004	return find1 (0, 0, 0);
1005	}
1006	BOOST_UBLAS_INLINE
1007	iterator1 end1 () {
1008	return find1 (0, size1_, 0);
1009	}
1010
1011	class const_iterator2:
1012	public container_const_reference<mapped_matrix>,
1013	public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
1014	const_iterator2, value_type> {
1015	public:
1016	typedef typename mapped_matrix::value_type value_type;
1017	typedef typename mapped_matrix::difference_type difference_type;
1018	typedef typename mapped_matrix::const_reference reference;
1019	typedef const typename mapped_matrix::pointer pointer;
1020
1021	typedef const_iterator1 dual_iterator_type;
1022	typedef const_reverse_iterator1 dual_reverse_iterator_type;
1023
1024	// Construction and destruction
1025	BOOST_UBLAS_INLINE
1026	const_iterator2 ():
1027	container_const_reference<self_type> (), rank_ (), i_ (), j_ (), it_ () {}
1028	BOOST_UBLAS_INLINE
1029	const_iterator2 (const self_type &m, int rank, size_type i, size_type j, const const_subiterator_type &it):
1030	container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), it_ (it) {}
1031	BOOST_UBLAS_INLINE
1032	const_iterator2 (const iterator2 &it):
1033	container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), it_ (it.it_) {}
1034
1035	// Arithmetic
1036	BOOST_UBLAS_INLINE
1037	const_iterator2 &operator ++ () {
1038	if (rank_ == 1 && layout_type::fast_j ())
1039	++ it_;
1040	else
1041	*this = (*this) ().find2 (rank_, i_, index2 () + 1, 1);
1042	return *this;
1043	}
1044	BOOST_UBLAS_INLINE
1045	const_iterator2 &operator -- () {
1046	if (rank_ == 1 && layout_type::fast_j ())
1047	-- it_;
1048	else
1049	*this = (*this) ().find2 (rank_, i_, index2 () - 1, -1);
1050	return *this;
1051	}
1052
1053	// Dereference
1054	BOOST_UBLAS_INLINE
1055	const_reference operator * () const {
1056	BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
1057	BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
1058	if (rank_ == 1) {
1059	return (*it_).second;
1060	} else {
1061	return (*this) () (i_, j_);
1062	}
1063	}
1064
1065	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1066	BOOST_UBLAS_INLINE
1067	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1068	typename self_type::
1069	#endif
1070	const_iterator1 begin () const {
1071	const self_type &m = (*this) ();
1072	return m.find1 (1, 0, index2 ());
1073	}
1074	BOOST_UBLAS_INLINE
1075	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1076	typename self_type::
1077	#endif
1078	const_iterator1 cbegin () const {
1079	return begin ();
1080	}
1081	BOOST_UBLAS_INLINE
1082	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1083	typename self_type::
1084	#endif
1085	const_iterator1 end () const {
1086	const self_type &m = (*this) ();
1087	return m.find1 (1, m.size1 (), index2 ());
1088	}
1089	BOOST_UBLAS_INLINE
1090	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1091	typename self_type::
1092	#endif
1093	const_iterator1 cend () const {
1094	return end ();
1095	}
1096	BOOST_UBLAS_INLINE
1097	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1098	typename self_type::
1099	#endif
1100	const_reverse_iterator1 rbegin () const {
1101	return const_reverse_iterator1 (end ());
1102	}
1103	BOOST_UBLAS_INLINE
1104	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1105	typename self_type::
1106	#endif
1107	const_reverse_iterator1 crbegin () const {
1108	return rbegin ();
1109	}
1110	BOOST_UBLAS_INLINE
1111	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1112	typename self_type::
1113	#endif
1114	const_reverse_iterator1 rend () const {
1115	return const_reverse_iterator1 (begin ());
1116	}
1117	BOOST_UBLAS_INLINE
1118	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1119	typename self_type::
1120	#endif
1121	const_reverse_iterator1 crend () const {
1122	return rend ();
1123	}
1124	#endif
1125
1126	// Indices
1127	BOOST_UBLAS_INLINE
1128	size_type index1 () const {
1129	if (rank_ == 1) {
1130	const self_type &m = (*this) ();
1131	BOOST_UBLAS_CHECK (layout_type::index_i ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size1 (), bad_index ());
1132	return layout_type::index_i ((*it_).first, m.size1 (), m.size2 ());
1133	} else {
1134	return i_;
1135	}
1136	}
1137	BOOST_UBLAS_INLINE
1138	size_type index2 () const {
1139	BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
1140	if (rank_ == 1) {
1141	const self_type &m = (*this) ();
1142	BOOST_UBLAS_CHECK (layout_type::index_j ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size2 (), bad_index ());
1143	return layout_type::index_j ((*it_).first, m.size1 (), m.size2 ());
1144	} else {
1145	return j_;
1146	}
1147	}
1148
1149	// Assignment
1150	BOOST_UBLAS_INLINE
1151	const_iterator2 &operator = (const const_iterator2 &it) {
1152	container_const_reference<self_type>::assign (&it ());
1153	rank_ = it.rank_;
1154	i_ = it.i_;
1155	j_ = it.j_;
1156	it_ = it.it_;
1157	return *this;
1158	}
1159
1160	// Comparison
1161	BOOST_UBLAS_INLINE
1162	bool operator == (const const_iterator2 &it) const {
1163	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1164	// BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
1165	if (rank_ == 1 || it.rank_ == 1) {
1166	return it_ == it.it_;
1167	} else {
1168	return i_ == it.i_ && j_ == it.j_;
1169	}
1170	}
1171
1172	private:
1173	int rank_;
1174	size_type i_;
1175	size_type j_;
1176	const_subiterator_type it_;
1177	};
1178
1179	BOOST_UBLAS_INLINE
1180	const_iterator2 begin2 () const {
1181	return find2 (0, 0, 0);
1182	}
1183	BOOST_UBLAS_INLINE
1184	const_iterator2 cbegin2 () const {
1185	return begin2 ();
1186	}
1187	BOOST_UBLAS_INLINE
1188	const_iterator2 end2 () const {
1189	return find2 (0, 0, size2_);
1190	}
1191	BOOST_UBLAS_INLINE
1192	const_iterator2 cend2 () const {
1193	return end2 ();
1194	}
1195
1196	class iterator2:
1197	public container_reference<mapped_matrix>,
1198	public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
1199	iterator2, value_type> {
1200	public:
1201	typedef typename mapped_matrix::value_type value_type;
1202	typedef typename mapped_matrix::difference_type difference_type;
1203	typedef typename mapped_matrix::true_reference reference;
1204	typedef typename mapped_matrix::pointer pointer;
1205
1206	typedef iterator1 dual_iterator_type;
1207	typedef reverse_iterator1 dual_reverse_iterator_type;
1208
1209	// Construction and destruction
1210	BOOST_UBLAS_INLINE
1211	iterator2 ():
1212	container_reference<self_type> (), rank_ (), i_ (), j_ (), it_ () {}
1213	BOOST_UBLAS_INLINE
1214	iterator2 (self_type &m, int rank, size_type i, size_type j, const subiterator_type &it):
1215	container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), it_ (it) {}
1216
1217	// Arithmetic
1218	BOOST_UBLAS_INLINE
1219	iterator2 &operator ++ () {
1220	if (rank_ == 1 && layout_type::fast_j ())
1221	++ it_;
1222	else
1223	*this = (*this) ().find2 (rank_, i_, index2 () + 1, 1);
1224	return *this;
1225	}
1226	BOOST_UBLAS_INLINE
1227	iterator2 &operator -- () {
1228	if (rank_ == 1 && layout_type::fast_j ())
1229	-- it_;
1230	else
1231	*this = (*this) ().find2 (rank_, i_, index2 () - 1, -1);
1232	return *this;
1233	}
1234
1235	// Dereference
1236	BOOST_UBLAS_INLINE
1237	reference operator * () const {
1238	BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
1239	BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
1240	if (rank_ == 1) {
1241	return (*it_).second;
1242	} else {
1243	return (*this) ().at_element (i_, j_);
1244	}
1245	}
1246
1247	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1248	BOOST_UBLAS_INLINE
1249	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1250	typename self_type::
1251	#endif
1252	iterator1 begin () const {
1253	self_type &m = (*this) ();
1254	return m.find1 (1, 0, index2 ());
1255	}
1256	BOOST_UBLAS_INLINE
1257	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1258	typename self_type::
1259	#endif
1260	iterator1 end () const {
1261	self_type &m = (*this) ();
1262	return m.find1 (1, m.size1 (), index2 ());
1263	}
1264	BOOST_UBLAS_INLINE
1265	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1266	typename self_type::
1267	#endif
1268	reverse_iterator1 rbegin () const {
1269	return reverse_iterator1 (end ());
1270	}
1271	BOOST_UBLAS_INLINE
1272	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1273	typename self_type::
1274	#endif
1275	reverse_iterator1 rend () const {
1276	return reverse_iterator1 (begin ());
1277	}
1278	#endif
1279
1280	// Indices
1281	BOOST_UBLAS_INLINE
1282	size_type index1 () const {
1283	if (rank_ == 1) {
1284	const self_type &m = (*this) ();
1285	BOOST_UBLAS_CHECK (layout_type::index_i ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size1 (), bad_index ());
1286	return layout_type::index_i ((*it_).first, m.size1 (), m.size2 ());
1287	} else {
1288	return i_;
1289	}
1290	}
1291	BOOST_UBLAS_INLINE
1292	size_type index2 () const {
1293	BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
1294	if (rank_ == 1) {
1295	const self_type &m = (*this) ();
1296	BOOST_UBLAS_CHECK (layout_type::index_j ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size2 (), bad_index ());
1297	return layout_type::index_j ((*it_).first, m.size1 (), m.size2 ());
1298	} else {
1299	return j_;
1300	}
1301	}
1302
1303	// Assignment
1304	BOOST_UBLAS_INLINE
1305	iterator2 &operator = (const iterator2 &it) {
1306	container_reference<self_type>::assign (&it ());
1307	rank_ = it.rank_;
1308	i_ = it.i_;
1309	j_ = it.j_;
1310	it_ = it.it_;
1311	return *this;
1312	}
1313
1314	// Comparison
1315	BOOST_UBLAS_INLINE
1316	bool operator == (const iterator2 &it) const {
1317	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1318	// BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
1319	if (rank_ == 1 || it.rank_ == 1) {
1320	return it_ == it.it_;
1321	} else {
1322	return i_ == it.i_ && j_ == it.j_;
1323	}
1324	}
1325
1326	private:
1327	int rank_;
1328	size_type i_;
1329	size_type j_;
1330	subiterator_type it_;
1331
1332	friend class const_iterator2;
1333	};
1334
1335	BOOST_UBLAS_INLINE
1336	iterator2 begin2 () {
1337	return find2 (0, 0, 0);
1338	}
1339	BOOST_UBLAS_INLINE
1340	iterator2 end2 () {
1341	return find2 (0, 0, size2_);
1342	}
1343
1344	// Reverse iterators
1345
1346	BOOST_UBLAS_INLINE
1347	const_reverse_iterator1 rbegin1 () const {
1348	return const_reverse_iterator1 (end1 ());
1349	}
1350	BOOST_UBLAS_INLINE
1351	const_reverse_iterator1 crbegin1 () const {
1352	return rbegin1 ();
1353	}
1354	BOOST_UBLAS_INLINE
1355	const_reverse_iterator1 rend1 () const {
1356	return const_reverse_iterator1 (begin1 ());
1357	}
1358	BOOST_UBLAS_INLINE
1359	const_reverse_iterator1 crend1 () const {
1360	return rend1 ();
1361	}
1362
1363	BOOST_UBLAS_INLINE
1364	reverse_iterator1 rbegin1 () {
1365	return reverse_iterator1 (end1 ());
1366	}
1367	BOOST_UBLAS_INLINE
1368	reverse_iterator1 rend1 () {
1369	return reverse_iterator1 (begin1 ());
1370	}
1371
1372	BOOST_UBLAS_INLINE
1373	const_reverse_iterator2 rbegin2 () const {
1374	return const_reverse_iterator2 (end2 ());
1375	}
1376	BOOST_UBLAS_INLINE
1377	const_reverse_iterator2 crbegin2 () const {
1378	return rbegin2 ();
1379	}
1380	BOOST_UBLAS_INLINE
1381	const_reverse_iterator2 rend2 () const {
1382	return const_reverse_iterator2 (begin2 ());
1383	}
1384	BOOST_UBLAS_INLINE
1385	const_reverse_iterator2 crend2 () const {
1386	return rend2 ();
1387	}
1388
1389	BOOST_UBLAS_INLINE
1390	reverse_iterator2 rbegin2 () {
1391	return reverse_iterator2 (end2 ());
1392	}
1393	BOOST_UBLAS_INLINE
1394	reverse_iterator2 rend2 () {
1395	return reverse_iterator2 (begin2 ());
1396	}
1397
1398	// Serialization
1399	template<class Archive>
1400	void serialize(Archive & ar, const unsigned int /* file_version */){
1401	serialization::collection_size_type s1 (size1_);
1402	serialization::collection_size_type s2 (size2_);
1403	ar & serialization::make_nvp(n: "size1",v&: s1);
1404	ar & serialization::make_nvp(n: "size2",v&: s2);
1405	if (Archive::is_loading::value) {
1406	size1_ = s1;
1407	size2_ = s2;
1408	}
1409	ar & serialization::make_nvp("data", data_);
1410	}
1411
1412	private:
1413	size_type size1_;
1414	size_type size2_;
1415	array_type data_;
1416	static const value_type zero_;
1417	};
1418
1419	template<class T, class L, class A>
1420	const typename mapped_matrix<T, L, A>::value_type mapped_matrix<T, L, A>::zero_ = value_type/*zero*/();
1421
1422
1423	// Vector index map based sparse matrix class
1424	template<class T, class L, class A>
1425	class mapped_vector_of_mapped_vector:
1426	public matrix_container<mapped_vector_of_mapped_vector<T, L, A> > {
1427
1428	typedef T &true_reference;
1429	typedef T *pointer;
1430	typedef const T *const_pointer;
1431	typedef A array_type;
1432	typedef const A const_array_type;
1433	typedef L layout_type;
1434	typedef mapped_vector_of_mapped_vector<T, L, A> self_type;
1435	public:
1436	#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
1437	using matrix_container<self_type>::operator ();
1438	#endif
1439	typedef typename A::size_type size_type;
1440	typedef typename A::difference_type difference_type;
1441	typedef T value_type;
1442	typedef const T &const_reference;
1443	#ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
1444	typedef typename detail::map_traits<typename A::data_value_type, T>::reference reference;
1445	#else
1446	typedef sparse_matrix_element<self_type> reference;
1447	#endif
1448	typedef const matrix_reference<const self_type> const_closure_type;
1449	typedef matrix_reference<self_type> closure_type;
1450	typedef mapped_vector<T> vector_temporary_type;
1451	typedef self_type matrix_temporary_type;
1452	typedef typename A::value_type::second_type vector_data_value_type;
1453	typedef sparse_tag storage_category;
1454	typedef typename L::orientation_category orientation_category;
1455
1456	// Construction and destruction
1457	BOOST_UBLAS_INLINE
1458	mapped_vector_of_mapped_vector ():
1459	matrix_container<self_type> (),
1460	size1_ (0), size2_ (0), data_ () {
1461	data_ [layout_type::size_M (size1_, size2_)] = vector_data_value_type ();
1462	}
1463	BOOST_UBLAS_INLINE
1464	mapped_vector_of_mapped_vector (size_type size1, size_type size2, size_type /*non_zeros*/ = 0):
1465	matrix_container<self_type> (),
1466	size1_ (size1), size2_ (size2), data_ () {
1467	data_ [layout_type::size_M (size1_, size2_)] = vector_data_value_type ();
1468	}
1469	BOOST_UBLAS_INLINE
1470	mapped_vector_of_mapped_vector (const mapped_vector_of_mapped_vector &m):
1471	matrix_container<self_type> (),
1472	size1_ (m.size1_), size2_ (m.size2_), data_ (m.data_) {}
1473	template<class AE>
1474	BOOST_UBLAS_INLINE
1475	mapped_vector_of_mapped_vector (const matrix_expression<AE> &ae, size_type /*non_zeros*/ = 0):
1476	matrix_container<self_type> (),
1477	size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), data_ () {
1478	data_ [layout_type::size_M (size1_, size2_)] = vector_data_value_type ();
1479	matrix_assign<scalar_assign> (*this, ae);
1480	}
1481
1482	// Accessors
1483	BOOST_UBLAS_INLINE
1484	size_type size1 () const {
1485	return size1_;
1486	}
1487	BOOST_UBLAS_INLINE
1488	size_type size2 () const {
1489	return size2_;
1490	}
1491	BOOST_UBLAS_INLINE
1492	size_type nnz_capacity () const {
1493	size_type non_zeros = 0;
1494	for (vector_const_subiterator_type itv = data_ ().begin (); itv != data_ ().end (); ++ itv)
1495	non_zeros += detail::map_capacity (*itv);
1496	return non_zeros;
1497	}
1498	BOOST_UBLAS_INLINE
1499	size_type nnz () const {
1500	size_type filled = 0;
1501	for (vector_const_subiterator_type itv = data_ ().begin (); itv != data_ ().end (); ++ itv)
1502	filled += (*itv).size ();
1503	return filled;
1504	}
1505
1506	// Storage accessors
1507	BOOST_UBLAS_INLINE
1508	const_array_type &data () const {
1509	return data_;
1510	}
1511	BOOST_UBLAS_INLINE
1512	array_type &data () {
1513	return data_;
1514	}
1515
1516	// Resizing
1517	BOOST_UBLAS_INLINE
1518	void resize (size_type size1, size_type size2, bool preserve = true) {
1519	// FIXME preserve unimplemented
1520	BOOST_UBLAS_CHECK (!preserve, internal_logic ());
1521	size1_ = size1;
1522	size2_ = size2;
1523	data ().clear ();
1524	data () [layout_type::size_M (size1_, size2_)] = vector_data_value_type ();
1525	}
1526
1527	// Element support
1528	BOOST_UBLAS_INLINE
1529	pointer find_element (size_type i, size_type j) {
1530	return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i, j));
1531	}
1532	BOOST_UBLAS_INLINE
1533	const_pointer find_element (size_type i, size_type j) const {
1534	const size_type element1 = layout_type::index_M (i, j);
1535	const size_type element2 = layout_type::index_m (i, j);
1536	vector_const_subiterator_type itv (data ().find (element1));
1537	if (itv == data ().end ())
1538	return 0;
1539	BOOST_UBLAS_CHECK ((*itv).first == element1, internal_logic ()); // broken map
1540	const_subiterator_type it ((*itv).second.find (element2));
1541	if (it == (*itv).second.end ())
1542	return 0;
1543	BOOST_UBLAS_CHECK ((*it).first == element2, internal_logic ()); // broken map
1544	return &(*it).second;
1545	}
1546
1547	// Element access
1548	BOOST_UBLAS_INLINE
1549	const_reference operator () (size_type i, size_type j) const {
1550	const size_type element1 = layout_type::index_M (i, j);
1551	const size_type element2 = layout_type::index_m (i, j);
1552	vector_const_subiterator_type itv (data ().find (element1));
1553	if (itv == data ().end ())
1554	return zero_;
1555	BOOST_UBLAS_CHECK ((*itv).first == element1, internal_logic ()); // broken map
1556	const_subiterator_type it ((*itv).second.find (element2));
1557	if (it == (*itv).second.end ())
1558	return zero_;
1559	BOOST_UBLAS_CHECK ((*itv).first == element1, internal_logic ()); // broken map
1560	return (*it).second;
1561	}
1562	BOOST_UBLAS_INLINE
1563	reference operator () (size_type i, size_type j) {
1564	#ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
1565	const size_type element1 = layout_type::index_M (i, j);
1566	const size_type element2 = layout_type::index_m (i, j);
1567	vector_data_value_type& vd (data () [element1]);
1568	std::pair<subiterator_type, bool> ii (vd.insert (typename array_type::value_type::second_type::value_type (element2, value_type/*zero*/())));
1569	BOOST_UBLAS_CHECK ((ii.first)->first == element2, internal_logic ()); // broken map
1570	return (ii.first)->second;
1571	#else
1572	return reference (*this, i, j);
1573	#endif
1574	}
1575
1576	// Element assignment
1577	BOOST_UBLAS_INLINE
1578	true_reference insert_element (size_type i, size_type j, const_reference t) {
1579	BOOST_UBLAS_CHECK (!find_element (i, j), bad_index ()); // duplicate element
1580	const size_type element1 = layout_type::index_M (i, j);
1581	const size_type element2 = layout_type::index_m (i, j);
1582
1583	vector_data_value_type& vd (data () [element1]);
1584	std::pair<subiterator_type, bool> ii (vd.insert (typename vector_data_value_type::value_type (element2, t)));
1585	BOOST_UBLAS_CHECK ((ii.first)->first == element2, internal_logic ()); // broken map
1586	if (!ii.second) // existing element
1587	(ii.first)->second = t;
1588	return (ii.first)->second;
1589	}
1590	BOOST_UBLAS_INLINE
1591	void erase_element (size_type i, size_type j) {
1592	vector_subiterator_type itv (data ().find (layout_type::index_M (i, j)));
1593	if (itv == data ().end ())
1594	return;
1595	subiterator_type it ((*itv).second.find (layout_type::index_m (i, j)));
1596	if (it == (*itv).second.end ())
1597	return;
1598	(*itv).second.erase (it);
1599	}
1600
1601	// Zeroing
1602	BOOST_UBLAS_INLINE
1603	void clear () {
1604	data ().clear ();
1605	data_ [layout_type::size_M (size1_, size2_)] = vector_data_value_type ();
1606	}
1607
1608	// Assignment
1609	BOOST_UBLAS_INLINE
1610	mapped_vector_of_mapped_vector &operator = (const mapped_vector_of_mapped_vector &m) {
1611	if (this != &m) {
1612	size1_ = m.size1_;
1613	size2_ = m.size2_;
1614	data () = m.data ();
1615	}
1616	return *this;
1617	}
1618	template<class C> // Container assignment without temporary
1619	BOOST_UBLAS_INLINE
1620	mapped_vector_of_mapped_vector &operator = (const matrix_container<C> &m) {
1621	resize (size1: m ().size1 (), size2: m ().size2 (), preserve: false);
1622	assign (m);
1623	return *this;
1624	}
1625	BOOST_UBLAS_INLINE
1626	mapped_vector_of_mapped_vector &assign_temporary (mapped_vector_of_mapped_vector &m) {
1627	swap (m);
1628	return *this;
1629	}
1630	template<class AE>
1631	BOOST_UBLAS_INLINE
1632	mapped_vector_of_mapped_vector &operator = (const matrix_expression<AE> &ae) {
1633	self_type temporary (ae);
1634	return assign_temporary (m&: temporary);
1635	}
1636	template<class AE>
1637	BOOST_UBLAS_INLINE
1638	mapped_vector_of_mapped_vector &assign (const matrix_expression<AE> &ae) {
1639	matrix_assign<scalar_assign> (*this, ae);
1640	return *this;
1641	}
1642	template<class AE>
1643	BOOST_UBLAS_INLINE
1644	mapped_vector_of_mapped_vector& operator += (const matrix_expression<AE> &ae) {
1645	self_type temporary (*this + ae);
1646	return assign_temporary (m&: temporary);
1647	}
1648	template<class C> // Container assignment without temporary
1649	BOOST_UBLAS_INLINE
1650	mapped_vector_of_mapped_vector &operator += (const matrix_container<C> &m) {
1651	plus_assign (m);
1652	return *this;
1653	}
1654	template<class AE>
1655	BOOST_UBLAS_INLINE
1656	mapped_vector_of_mapped_vector &plus_assign (const matrix_expression<AE> &ae) {
1657	matrix_assign<scalar_plus_assign> (*this, ae);
1658	return *this;
1659	}
1660	template<class AE>
1661	BOOST_UBLAS_INLINE
1662	mapped_vector_of_mapped_vector& operator -= (const matrix_expression<AE> &ae) {
1663	self_type temporary (*this - ae);
1664	return assign_temporary (m&: temporary);
1665	}
1666	template<class C> // Container assignment without temporary
1667	BOOST_UBLAS_INLINE
1668	mapped_vector_of_mapped_vector &operator -= (const matrix_container<C> &m) {
1669	minus_assign (m);
1670	return *this;
1671	}
1672	template<class AE>
1673	BOOST_UBLAS_INLINE
1674	mapped_vector_of_mapped_vector &minus_assign (const matrix_expression<AE> &ae) {
1675	matrix_assign<scalar_minus_assign> (*this, ae);
1676	return *this;
1677	}
1678	template<class AT>
1679	BOOST_UBLAS_INLINE
1680	mapped_vector_of_mapped_vector& operator *= (const AT &at) {
1681	matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
1682	return *this;
1683	}
1684	template<class AT>
1685	BOOST_UBLAS_INLINE
1686	mapped_vector_of_mapped_vector& operator /= (const AT &at) {
1687	matrix_assign_scalar<scalar_divides_assign> (*this, at);
1688	return *this;
1689	}
1690
1691	// Swapping
1692	BOOST_UBLAS_INLINE
1693	void swap (mapped_vector_of_mapped_vector &m) {
1694	if (this != &m) {
1695	std::swap (size1_, m.size1_);
1696	std::swap (size2_, m.size2_);
1697	data ().swap (m.data ());
1698	}
1699	}
1700	BOOST_UBLAS_INLINE
1701	friend void swap (mapped_vector_of_mapped_vector &m1, mapped_vector_of_mapped_vector &m2) {
1702	m1.swap (m2);
1703	}
1704
1705	// Iterator types
1706	private:
1707	// Use storage iterators
1708	typedef typename A::const_iterator vector_const_subiterator_type;
1709	typedef typename A::iterator vector_subiterator_type;
1710	typedef typename A::value_type::second_type::const_iterator const_subiterator_type;
1711	typedef typename A::value_type::second_type::iterator subiterator_type;
1712
1713	BOOST_UBLAS_INLINE
1714	true_reference at_element (size_type i, size_type j) {
1715	const size_type element1 = layout_type::index_M (i, j);
1716	const size_type element2 = layout_type::index_m (i, j);
1717	vector_subiterator_type itv (data ().find (element1));
1718	BOOST_UBLAS_CHECK (itv != data ().end(), bad_index ());
1719	BOOST_UBLAS_CHECK ((*itv).first == element1, internal_logic ()); // broken map
1720	subiterator_type it ((*itv).second.find (element2));
1721	BOOST_UBLAS_CHECK (it != (*itv).second.end (), bad_index ());
1722	BOOST_UBLAS_CHECK ((*it).first == element2, internal_logic ()); // broken map
1723
1724	return it->second;
1725	}
1726
1727	public:
1728	class const_iterator1;
1729	class iterator1;
1730	class const_iterator2;
1731	class iterator2;
1732	typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
1733	typedef reverse_iterator_base1<iterator1> reverse_iterator1;
1734	typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
1735	typedef reverse_iterator_base2<iterator2> reverse_iterator2;
1736
1737	// Element lookup
1738	// BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1739	const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const {
1740	BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ());
1741	for (;;) {
1742	vector_const_subiterator_type itv (data ().lower_bound (layout_type::index_M (i, j)));
1743	vector_const_subiterator_type itv_end (data ().end ());
1744	if (itv == itv_end)
1745	return const_iterator1 (*this, rank, i, j, itv_end, (*(-- itv)).second.end ());
1746
1747	const_subiterator_type it ((*itv).second.lower_bound (layout_type::index_m (i, j)));
1748	const_subiterator_type it_end ((*itv).second.end ());
1749	if (rank == 0) {
1750	// advance to the first available major index
1751	size_type M = itv->first;
1752	size_type m;
1753	if (it != it_end) {
1754	m = it->first;
1755	} else {
1756	m = layout_type::size_m(size1_, size2_);
1757	}
1758	size_type first_i = layout_type::index_M(M,m);
1759	return const_iterator1 (*this, rank, first_i, j, itv, it);
1760	}
1761	if (it != it_end && (*it).first == layout_type::index_m (i, j))
1762	return const_iterator1 (*this, rank, i, j, itv, it);
1763	if (direction > 0) {
1764	if (layout_type::fast_i ()) {
1765	if (it == it_end)
1766	return const_iterator1 (*this, rank, i, j, itv, it);
1767	i = (*it).first;
1768	} else {
1769	if (i >= size1_)
1770	return const_iterator1 (*this, rank, i, j, itv, it);
1771	++ i;
1772	}
1773	} else /* if (direction < 0) */ {
1774	if (layout_type::fast_i ()) {
1775	if (it == (*itv).second.begin ())
1776	return const_iterator1 (*this, rank, i, j, itv, it);
1777	-- it;
1778	i = (*it).first;
1779	} else {
1780	if (i == 0)
1781	return const_iterator1 (*this, rank, i, j, itv, it);
1782	-- i;
1783	}
1784	}
1785	}
1786	}
1787	// BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1788	iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) {
1789	BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ());
1790	for (;;) {
1791	vector_subiterator_type itv (data ().lower_bound (layout_type::index_M (i, j)));
1792	vector_subiterator_type itv_end (data ().end ());
1793	if (itv == itv_end)
1794	return iterator1 (*this, rank, i, j, itv_end, (*(-- itv)).second.end ());
1795
1796	subiterator_type it ((*itv).second.lower_bound (layout_type::index_m (i, j)));
1797	subiterator_type it_end ((*itv).second.end ());
1798	if (rank == 0) {
1799	// advance to the first available major index
1800	size_type M = itv->first;
1801	size_type m;
1802	if (it != it_end) {
1803	m = it->first;
1804	} else {
1805	m = layout_type::size_m(size1_, size2_);
1806	}
1807	size_type first_i = layout_type::index_M(M,m);
1808	return iterator1 (*this, rank, first_i, j, itv, it);
1809	}
1810	if (it != it_end && (*it).first == layout_type::index_m (i, j))
1811	return iterator1 (*this, rank, i, j, itv, it);
1812	if (direction > 0) {
1813	if (layout_type::fast_i ()) {
1814	if (it == it_end)
1815	return iterator1 (*this, rank, i, j, itv, it);
1816	i = (*it).first;
1817	} else {
1818	if (i >= size1_)
1819	return iterator1 (*this, rank, i, j, itv, it);
1820	++ i;
1821	}
1822	} else /* if (direction < 0) */ {
1823	if (layout_type::fast_i ()) {
1824	if (it == (*itv).second.begin ())
1825	return iterator1 (*this, rank, i, j, itv, it);
1826	-- it;
1827	i = (*it).first;
1828	} else {
1829	if (i == 0)
1830	return iterator1 (*this, rank, i, j, itv, it);
1831	-- i;
1832	}
1833	}
1834	}
1835	}
1836	// BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1837	const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const {
1838	BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ());
1839	for (;;) {
1840	vector_const_subiterator_type itv (data ().lower_bound (layout_type::index_M (i, j)));
1841	vector_const_subiterator_type itv_end (data ().end ());
1842	if (itv == itv_end)
1843	return const_iterator2 (*this, rank, i, j, itv_end, (*(-- itv)).second.end ());
1844
1845	const_subiterator_type it ((*itv).second.lower_bound (layout_type::index_m (i, j)));
1846	const_subiterator_type it_end ((*itv).second.end ());
1847	if (rank == 0) {
1848	// advance to the first available major index
1849	size_type M = itv->first;
1850	size_type m;
1851	if (it != it_end) {
1852	m = it->first;
1853	} else {
1854	m = layout_type::size_m(size1_, size2_);
1855	}
1856	size_type first_j = layout_type::index_m(M,m);
1857	return const_iterator2 (*this, rank, i, first_j, itv, it);
1858	}
1859	if (it != it_end && (*it).first == layout_type::index_m (i, j))
1860	return const_iterator2 (*this, rank, i, j, itv, it);
1861	if (direction > 0) {
1862	if (layout_type::fast_j ()) {
1863	if (it == it_end)
1864	return const_iterator2 (*this, rank, i, j, itv, it);
1865	j = (*it).first;
1866	} else {
1867	if (j >= size2_)
1868	return const_iterator2 (*this, rank, i, j, itv, it);
1869	++ j;
1870	}
1871	} else /* if (direction < 0) */ {
1872	if (layout_type::fast_j ()) {
1873	if (it == (*itv).second.begin ())
1874	return const_iterator2 (*this, rank, i, j, itv, it);
1875	-- it;
1876	j = (*it).first;
1877	} else {
1878	if (j == 0)
1879	return const_iterator2 (*this, rank, i, j, itv, it);
1880	-- j;
1881	}
1882	}
1883	}
1884	}
1885	// BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1886	iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) {
1887	BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ());
1888	for (;;) {
1889	vector_subiterator_type itv (data ().lower_bound (layout_type::index_M (i, j)));
1890	vector_subiterator_type itv_end (data ().end ());
1891	if (itv == itv_end)
1892	return iterator2 (*this, rank, i, j, itv_end, (*(-- itv)).second.end ());
1893
1894	subiterator_type it ((*itv).second.lower_bound (layout_type::index_m (i, j)));
1895	subiterator_type it_end ((*itv).second.end ());
1896	if (rank == 0) {
1897	// advance to the first available major index
1898	size_type M = itv->first;
1899	size_type m;
1900	if (it != it_end) {
1901	m = it->first;
1902	} else {
1903	m = layout_type::size_m(size1_, size2_);
1904	}
1905	size_type first_j = layout_type::index_m(M,m);
1906	return iterator2 (*this, rank, i, first_j, itv, it);
1907	}
1908	if (it != it_end && (*it).first == layout_type::index_m (i, j))
1909	return iterator2 (*this, rank, i, j, itv, it);
1910	if (direction > 0) {
1911	if (layout_type::fast_j ()) {
1912	if (it == it_end)
1913	return iterator2 (*this, rank, i, j, itv, it);
1914	j = (*it).first;
1915	} else {
1916	if (j >= size2_)
1917	return iterator2 (*this, rank, i, j, itv, it);
1918	++ j;
1919	}
1920	} else /* if (direction < 0) */ {
1921	if (layout_type::fast_j ()) {
1922	if (it == (*itv).second.begin ())
1923	return iterator2 (*this, rank, i, j, itv, it);
1924	-- it;
1925	j = (*it).first;
1926	} else {
1927	if (j == 0)
1928	return iterator2 (*this, rank, i, j, itv, it);
1929	-- j;
1930	}
1931	}
1932	}
1933	}
1934
1935	class const_iterator1:
1936	public container_const_reference<mapped_vector_of_mapped_vector>,
1937	public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
1938	const_iterator1, value_type> {
1939	public:
1940	typedef typename mapped_vector_of_mapped_vector::value_type value_type;
1941	typedef typename mapped_vector_of_mapped_vector::difference_type difference_type;
1942	typedef typename mapped_vector_of_mapped_vector::const_reference reference;
1943	typedef const typename mapped_vector_of_mapped_vector::pointer pointer;
1944
1945	typedef const_iterator2 dual_iterator_type;
1946	typedef const_reverse_iterator2 dual_reverse_iterator_type;
1947
1948	// Construction and destruction
1949	BOOST_UBLAS_INLINE
1950	const_iterator1 ():
1951	container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
1952	BOOST_UBLAS_INLINE
1953	const_iterator1 (const self_type &m, int rank, size_type i, size_type j, const vector_const_subiterator_type &itv, const const_subiterator_type &it):
1954	container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
1955	BOOST_UBLAS_INLINE
1956	const_iterator1 (const iterator1 &it):
1957	container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {}
1958
1959	// Arithmetic
1960	BOOST_UBLAS_INLINE
1961	const_iterator1 &operator ++ () {
1962	if (rank_ == 1 && layout_type::fast_i ())
1963	++ it_;
1964	else {
1965	const self_type &m = (*this) ();
1966	if (rank_ == 0) {
1967	++ itv_;
1968	i_ = itv_->first;
1969	} else {
1970	i_ = index1 () + 1;
1971	}
1972	if (rank_ == 1 && ++ itv_ == m.end1 ().itv_)
1973	*this = m.find1 (rank_, i_, j_, 1);
1974	else if (rank_ == 1) {
1975	it_ = (*itv_).second.begin ();
1976	if (it_ == (*itv_).second.end () || index2 () != j_)
1977	*this = m.find1 (rank_, i_, j_, 1);
1978	}
1979	}
1980	return *this;
1981	}
1982	BOOST_UBLAS_INLINE
1983	const_iterator1 &operator -- () {
1984	if (rank_ == 1 && layout_type::fast_i ())
1985	-- it_;
1986	else {
1987	const self_type &m = (*this) ();
1988	if (rank_ == 0) {
1989	-- itv_;
1990	i_ = itv_->first;
1991	} else {
1992	i_ = index1 () - 1;
1993	}
1994	// FIXME: this expression should never become true!
1995	if (rank_ == 1 && -- itv_ == m.end1 ().itv_)
1996	*this = m.find1 (rank_, i_, j_, -1);
1997	else if (rank_ == 1) {
1998	it_ = (*itv_).second.begin ();
1999	if (it_ == (*itv_).second.end () || index2 () != j_)
2000	*this = m.find1 (rank_, i_, j_, -1);
2001	}
2002	}
2003	return *this;
2004	}
2005
2006	// Dereference
2007	BOOST_UBLAS_INLINE
2008	const_reference operator * () const {
2009	BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
2010	BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
2011	if (rank_ == 1) {
2012	return (*it_).second;
2013	} else {
2014	return (*this) () (i_, j_);
2015	}
2016	}
2017
2018	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
2019	BOOST_UBLAS_INLINE
2020	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2021	typename self_type::
2022	#endif
2023	const_iterator2 begin () const {
2024	const self_type &m = (*this) ();
2025	return m.find2 (1, index1 (), 0);
2026	}
2027	BOOST_UBLAS_INLINE
2028	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2029	typename self_type::
2030	#endif
2031	const_iterator2 cbegin () const {
2032	return begin ();
2033	}
2034	BOOST_UBLAS_INLINE
2035	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2036	typename self_type::
2037	#endif
2038	const_iterator2 end () const {
2039	const self_type &m = (*this) ();
2040	return m.find2 (1, index1 (), m.size2 ());
2041	}
2042	BOOST_UBLAS_INLINE
2043	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2044	typename self_type::
2045	#endif
2046	const_iterator2 cend () const {
2047	return end ();
2048	}
2049	BOOST_UBLAS_INLINE
2050	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2051	typename self_type::
2052	#endif
2053	const_reverse_iterator2 rbegin () const {
2054	return const_reverse_iterator2 (end ());
2055	}
2056	BOOST_UBLAS_INLINE
2057	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2058	typename self_type::
2059	#endif
2060	const_reverse_iterator2 crbegin () const {
2061	return rbegin ();
2062	}
2063	BOOST_UBLAS_INLINE
2064	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2065	typename self_type::
2066	#endif
2067	const_reverse_iterator2 rend () const {
2068	return const_reverse_iterator2 (begin ());
2069	}
2070	BOOST_UBLAS_INLINE
2071	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2072	typename self_type::
2073	#endif
2074	const_reverse_iterator2 crend () const {
2075	return rend ();
2076	}
2077	#endif
2078
2079	// Indices
2080	BOOST_UBLAS_INLINE
2081	size_type index1 () const {
2082	BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
2083	if (rank_ == 1) {
2084	BOOST_UBLAS_CHECK (layout_type::index_M ((*itv_).first, (*it_).first) < (*this) ().size1 (), bad_index ());
2085	return layout_type::index_M ((*itv_).first, (*it_).first);
2086	} else {
2087	return i_;
2088	}
2089	}
2090	BOOST_UBLAS_INLINE
2091	size_type index2 () const {
2092	if (rank_ == 1) {
2093	BOOST_UBLAS_CHECK (layout_type::index_m ((*itv_).first, (*it_).first) < (*this) ().size2 (), bad_index ());
2094	return layout_type::index_m ((*itv_).first, (*it_).first);
2095	} else {
2096	return j_;
2097	}
2098	}
2099
2100	// Assignment
2101	BOOST_UBLAS_INLINE
2102	const_iterator1 &operator = (const const_iterator1 &it) {
2103	container_const_reference<self_type>::assign (&it ());
2104	rank_ = it.rank_;
2105	i_ = it.i_;
2106	j_ = it.j_;
2107	itv_ = it.itv_;
2108	it_ = it.it_;
2109	return *this;
2110	}
2111
2112	// Comparison
2113	BOOST_UBLAS_INLINE
2114	bool operator == (const const_iterator1 &it) const {
2115	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2116	// BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
2117	if (rank_ == 1 || it.rank_ == 1) {
2118	return it_ == it.it_;
2119	} else {
2120	return i_ == it.i_ && j_ == it.j_;
2121	}
2122	}
2123
2124	private:
2125	int rank_;
2126	size_type i_;
2127	size_type j_;
2128	vector_const_subiterator_type itv_;
2129	const_subiterator_type it_;
2130	};
2131
2132	BOOST_UBLAS_INLINE
2133	const_iterator1 begin1 () const {
2134	return find1 (0, 0, 0);
2135	}
2136	BOOST_UBLAS_INLINE
2137	const_iterator1 cbegin1 () const {
2138	return begin1 ();
2139	}
2140	BOOST_UBLAS_INLINE
2141	const_iterator1 end1 () const {
2142	return find1 (0, size1_, 0);
2143	}
2144	BOOST_UBLAS_INLINE
2145	const_iterator1 cend1 () const {
2146	return end1 ();
2147	}
2148
2149	class iterator1:
2150	public container_reference<mapped_vector_of_mapped_vector>,
2151	public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
2152	iterator1, value_type> {
2153	public:
2154	typedef typename mapped_vector_of_mapped_vector::value_type value_type;
2155	typedef typename mapped_vector_of_mapped_vector::difference_type difference_type;
2156	typedef typename mapped_vector_of_mapped_vector::true_reference reference;
2157	typedef typename mapped_vector_of_mapped_vector::pointer pointer;
2158
2159	typedef iterator2 dual_iterator_type;
2160	typedef reverse_iterator2 dual_reverse_iterator_type;
2161
2162	// Construction and destruction
2163	BOOST_UBLAS_INLINE
2164	iterator1 ():
2165	container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
2166	BOOST_UBLAS_INLINE
2167	iterator1 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it):
2168	container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
2169
2170	// Arithmetic
2171	BOOST_UBLAS_INLINE
2172	iterator1 &operator ++ () {
2173	if (rank_ == 1 && layout_type::fast_i ())
2174	++ it_;
2175	else {
2176	self_type &m = (*this) ();
2177	if (rank_ == 0) {
2178	++ itv_;
2179	i_ = itv_->first;
2180	} else {
2181	i_ = index1 () + 1;
2182	}
2183	if (rank_ == 1 && ++ itv_ == m.end1 ().itv_)
2184	*this = m.find1 (rank_, i_, j_, 1);
2185	else if (rank_ == 1) {
2186	it_ = (*itv_).second.begin ();
2187	if (it_ == (*itv_).second.end () || index2 () != j_)
2188	*this = m.find1 (rank_, i_, j_, 1);
2189	}
2190	}
2191	return *this;
2192	}
2193	BOOST_UBLAS_INLINE
2194	iterator1 &operator -- () {
2195	if (rank_ == 1 && layout_type::fast_i ())
2196	-- it_;
2197	else {
2198	self_type &m = (*this) ();
2199	if (rank_ == 0) {
2200	-- itv_;
2201	i_ = itv_->first;
2202	} else {
2203	i_ = index1 () - 1;
2204	}
2205	// FIXME: this expression should never become true!
2206	if (rank_ == 1 && -- itv_ == m.end1 ().itv_)
2207	*this = m.find1 (rank_, i_, j_, -1);
2208	else if (rank_ == 1) {
2209	it_ = (*itv_).second.begin ();
2210	if (it_ == (*itv_).second.end () || index2 () != j_)
2211	*this = m.find1 (rank_, i_, j_, -1);
2212	}
2213	}
2214	return *this;
2215	}
2216
2217	// Dereference
2218	BOOST_UBLAS_INLINE
2219	reference operator * () const {
2220	BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
2221	BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
2222	if (rank_ == 1) {
2223	return (*it_).second;
2224	} else {
2225	return (*this) ().at_element (i_, j_);
2226	}
2227	}
2228
2229	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
2230	BOOST_UBLAS_INLINE
2231	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2232	typename self_type::
2233	#endif
2234	iterator2 begin () const {
2235	self_type &m = (*this) ();
2236	return m.find2 (1, index1 (), 0);
2237	}
2238	BOOST_UBLAS_INLINE
2239	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2240	typename self_type::
2241	#endif
2242	iterator2 end () const {
2243	self_type &m = (*this) ();
2244	return m.find2 (1, index1 (), m.size2 ());
2245	}
2246	BOOST_UBLAS_INLINE
2247	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2248	typename self_type::
2249	#endif
2250	reverse_iterator2 rbegin () const {
2251	return reverse_iterator2 (end ());
2252	}
2253	BOOST_UBLAS_INLINE
2254	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2255	typename self_type::
2256	#endif
2257	reverse_iterator2 rend () const {
2258	return reverse_iterator2 (begin ());
2259	}
2260	#endif
2261
2262	// Indices
2263	BOOST_UBLAS_INLINE
2264	size_type index1 () const {
2265	BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
2266	if (rank_ == 1) {
2267	BOOST_UBLAS_CHECK (layout_type::index_M ((*itv_).first, (*it_).first) < (*this) ().size1 (), bad_index ());
2268	return layout_type::index_M ((*itv_).first, (*it_).first);
2269	} else {
2270	return i_;
2271	}
2272	}
2273	BOOST_UBLAS_INLINE
2274	size_type index2 () const {
2275	if (rank_ == 1) {
2276	BOOST_UBLAS_CHECK (layout_type::index_m ((*itv_).first, (*it_).first) < (*this) ().size2 (), bad_index ());
2277	return layout_type::index_m ((*itv_).first, (*it_).first);
2278	} else {
2279	return j_;
2280	}
2281	}
2282
2283	// Assignment
2284	BOOST_UBLAS_INLINE
2285	iterator1 &operator = (const iterator1 &it) {
2286	container_reference<self_type>::assign (&it ());
2287	rank_ = it.rank_;
2288	i_ = it.i_;
2289	j_ = it.j_;
2290	itv_ = it.itv_;
2291	it_ = it.it_;
2292	return *this;
2293	}
2294
2295	// Comparison
2296	BOOST_UBLAS_INLINE
2297	bool operator == (const iterator1 &it) const {
2298	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2299	// BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
2300	if (rank_ == 1 || it.rank_ == 1) {
2301	return it_ == it.it_;
2302	} else {
2303	return i_ == it.i_ && j_ == it.j_;
2304	}
2305	}
2306
2307	private:
2308	int rank_;
2309	size_type i_;
2310	size_type j_;
2311	vector_subiterator_type itv_;
2312	subiterator_type it_;
2313
2314	friend class const_iterator1;
2315	};
2316
2317	BOOST_UBLAS_INLINE
2318	iterator1 begin1 () {
2319	return find1 (0, 0, 0);
2320	}
2321	BOOST_UBLAS_INLINE
2322	iterator1 end1 () {
2323	return find1 (0, size1_, 0);
2324	}
2325
2326	class const_iterator2:
2327	public container_const_reference<mapped_vector_of_mapped_vector>,
2328	public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
2329	const_iterator2, value_type> {
2330	public:
2331	typedef typename mapped_vector_of_mapped_vector::value_type value_type;
2332	typedef typename mapped_vector_of_mapped_vector::difference_type difference_type;
2333	typedef typename mapped_vector_of_mapped_vector::const_reference reference;
2334	typedef const typename mapped_vector_of_mapped_vector::pointer pointer;
2335
2336	typedef const_iterator1 dual_iterator_type;
2337	typedef const_reverse_iterator1 dual_reverse_iterator_type;
2338
2339	// Construction and destruction
2340	BOOST_UBLAS_INLINE
2341	const_iterator2 ():
2342	container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
2343	BOOST_UBLAS_INLINE
2344	const_iterator2 (const self_type &m, int rank, size_type i, size_type j, const vector_const_subiterator_type &itv, const const_subiterator_type &it):
2345	container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
2346	BOOST_UBLAS_INLINE
2347	const_iterator2 (const iterator2 &it):
2348	container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {}
2349
2350	// Arithmetic
2351	BOOST_UBLAS_INLINE
2352	const_iterator2 &operator ++ () {
2353	if (rank_ == 1 && layout_type::fast_j ())
2354	++ it_;
2355	else {
2356	const self_type &m = (*this) ();
2357	if (rank_ == 0) {
2358	++ itv_;
2359	j_ = itv_->first;
2360	} else {
2361	j_ = index2 () + 1;
2362	}
2363	if (rank_ == 1 && ++ itv_ == m.end2 ().itv_)
2364	*this = m.find2 (rank_, i_, j_, 1);
2365	else if (rank_ == 1) {
2366	it_ = (*itv_).second.begin ();
2367	if (it_ == (*itv_).second.end () || index1 () != i_)
2368	*this = m.find2 (rank_, i_, j_, 1);
2369	}
2370	}
2371	return *this;
2372	}
2373	BOOST_UBLAS_INLINE
2374	const_iterator2 &operator -- () {
2375	if (rank_ == 1 && layout_type::fast_j ())
2376	-- it_;
2377	else {
2378	const self_type &m = (*this) ();
2379	if (rank_ == 0) {
2380	-- itv_;
2381	j_ = itv_->first;
2382	} else {
2383	j_ = index2 () - 1;
2384	}
2385	// FIXME: this expression should never become true!
2386	if (rank_ == 1 && -- itv_ == m.end2 ().itv_)
2387	*this = m.find2 (rank_, i_, j_, -1);
2388	else if (rank_ == 1) {
2389	it_ = (*itv_).second.begin ();
2390	if (it_ == (*itv_).second.end () || index1 () != i_)
2391	*this = m.find2 (rank_, i_, j_, -1);
2392	}
2393	}
2394	return *this;
2395	}
2396
2397	// Dereference
2398	BOOST_UBLAS_INLINE
2399	const_reference operator * () const {
2400	BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
2401	BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
2402	if (rank_ == 1) {
2403	return (*it_).second;
2404	} else {
2405	return (*this) () (i_, j_);
2406	}
2407	}
2408
2409	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
2410	BOOST_UBLAS_INLINE
2411	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2412	typename self_type::
2413	#endif
2414	const_iterator1 begin () const {
2415	const self_type &m = (*this) ();
2416	return m.find1 (1, 0, index2 ());
2417	}
2418	BOOST_UBLAS_INLINE
2419	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2420	typename self_type::
2421	#endif
2422	const_iterator1 cbegin () const {
2423	return begin ();
2424	}
2425	BOOST_UBLAS_INLINE
2426	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2427	typename self_type::
2428	#endif
2429	const_iterator1 end () const {
2430	const self_type &m = (*this) ();
2431	return m.find1 (1, m.size1 (), index2 ());
2432	}
2433	BOOST_UBLAS_INLINE
2434	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2435	typename self_type::
2436	#endif
2437	const_iterator1 cend () const {
2438	return end ();
2439	}
2440	BOOST_UBLAS_INLINE
2441	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2442	typename self_type::
2443	#endif
2444	const_reverse_iterator1 rbegin () const {
2445	return const_reverse_iterator1 (end ());
2446	}
2447	BOOST_UBLAS_INLINE
2448	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2449	typename self_type::
2450	#endif
2451	const_reverse_iterator1 crbegin () const {
2452	return rbegin ();
2453	}
2454	BOOST_UBLAS_INLINE
2455	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2456	typename self_type::
2457	#endif
2458	const_reverse_iterator1 rend () const {
2459	return const_reverse_iterator1 (begin ());
2460	}
2461	BOOST_UBLAS_INLINE
2462	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2463	typename self_type::
2464	#endif
2465	const_reverse_iterator1 crend () const {
2466	return rend ();
2467	}
2468	#endif
2469
2470	// Indices
2471	BOOST_UBLAS_INLINE
2472	size_type index1 () const {
2473	if (rank_ == 1) {
2474	BOOST_UBLAS_CHECK (layout_type::index_M ((*itv_).first, (*it_).first) < (*this) ().size1 (), bad_index ());
2475	return layout_type::index_M ((*itv_).first, (*it_).first);
2476	} else {
2477	return i_;
2478	}
2479	}
2480	BOOST_UBLAS_INLINE
2481	size_type index2 () const {
2482	BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
2483	if (rank_ == 1) {
2484	BOOST_UBLAS_CHECK (layout_type::index_m ((*itv_).first, (*it_).first) < (*this) ().size2 (), bad_index ());
2485	return layout_type::index_m ((*itv_).first, (*it_).first);
2486	} else {
2487	return j_;
2488	}
2489	}
2490
2491	// Assignment
2492	BOOST_UBLAS_INLINE
2493	const_iterator2 &operator = (const const_iterator2 &it) {
2494	container_const_reference<self_type>::assign (&it ());
2495	rank_ = it.rank_;
2496	i_ = it.i_;
2497	j_ = it.j_;
2498	itv_ = it.itv_;
2499	it_ = it.it_;
2500	return *this;
2501	}
2502
2503	// Comparison
2504	BOOST_UBLAS_INLINE
2505	bool operator == (const const_iterator2 &it) const {
2506	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2507	// BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
2508	if (rank_ == 1 || it.rank_ == 1) {
2509	return it_ == it.it_;
2510	} else {
2511	return i_ == it.i_ && j_ == it.j_;
2512	}
2513	}
2514
2515	private:
2516	int rank_;
2517	size_type i_;
2518	size_type j_;
2519	vector_const_subiterator_type itv_;
2520	const_subiterator_type it_;
2521	};
2522
2523	BOOST_UBLAS_INLINE
2524	const_iterator2 begin2 () const {
2525	return find2 (0, 0, 0);
2526	}
2527	BOOST_UBLAS_INLINE
2528	const_iterator2 cbegin2 () const {
2529	return begin2 ();
2530	}
2531	BOOST_UBLAS_INLINE
2532	const_iterator2 end2 () const {
2533	return find2 (0, 0, size2_);
2534	}
2535	BOOST_UBLAS_INLINE
2536	const_iterator2 cend2 () const {
2537	return end2 ();
2538	}
2539
2540	class iterator2:
2541	public container_reference<mapped_vector_of_mapped_vector>,
2542	public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
2543	iterator2, value_type> {
2544	public:
2545	typedef typename mapped_vector_of_mapped_vector::value_type value_type;
2546	typedef typename mapped_vector_of_mapped_vector::difference_type difference_type;
2547	typedef typename mapped_vector_of_mapped_vector::true_reference reference;
2548	typedef typename mapped_vector_of_mapped_vector::pointer pointer;
2549
2550	typedef iterator1 dual_iterator_type;
2551	typedef reverse_iterator1 dual_reverse_iterator_type;
2552
2553	// Construction and destruction
2554	BOOST_UBLAS_INLINE
2555	iterator2 ():
2556	container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
2557	BOOST_UBLAS_INLINE
2558	iterator2 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it):
2559	container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
2560
2561	// Arithmetic
2562	BOOST_UBLAS_INLINE
2563	iterator2 &operator ++ () {
2564	if (rank_ == 1 && layout_type::fast_j ())
2565	++ it_;
2566	else {
2567	self_type &m = (*this) ();
2568	if (rank_ == 0) {
2569	++ itv_;
2570	j_ = itv_->first;
2571	} else {
2572	j_ = index2 () + 1;
2573	}
2574	if (rank_ == 1 && ++ itv_ == m.end2 ().itv_)
2575	*this = m.find2 (rank_, i_, j_, 1);
2576	else if (rank_ == 1) {
2577	it_ = (*itv_).second.begin ();
2578	if (it_ == (*itv_).second.end () || index1 () != i_)
2579	*this = m.find2 (rank_, i_, j_, 1);
2580	}
2581	}
2582	return *this;
2583	}
2584	BOOST_UBLAS_INLINE
2585	iterator2 &operator -- () {
2586	if (rank_ == 1 && layout_type::fast_j ())
2587	-- it_;
2588	else {
2589	self_type &m = (*this) ();
2590	if (rank_ == 0) {
2591	-- itv_;
2592	j_ = itv_->first;
2593	} else {
2594	j_ = index2 () - 1;
2595	}
2596	// FIXME: this expression should never become true!
2597	if (rank_ == 1 && -- itv_ == m.end2 ().itv_)
2598	*this = m.find2 (rank_, i_, j_, -1);
2599	else if (rank_ == 1) {
2600	it_ = (*itv_).second.begin ();
2601	if (it_ == (*itv_).second.end () || index1 () != i_)
2602	*this = m.find2 (rank_, i_, j_, -1);
2603	}
2604	}
2605	return *this;
2606	}
2607
2608	// Dereference
2609	BOOST_UBLAS_INLINE
2610	reference operator * () const {
2611	BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
2612	BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
2613	if (rank_ == 1) {
2614	return (*it_).second;
2615	} else {
2616	return (*this) ().at_element (i_, j_);
2617	}
2618	}
2619
2620	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
2621	BOOST_UBLAS_INLINE
2622	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2623	typename self_type::
2624	#endif
2625	iterator1 begin () const {
2626	self_type &m = (*this) ();
2627	return m.find1 (1, 0, index2 ());
2628	}
2629	BOOST_UBLAS_INLINE
2630	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2631	typename self_type::
2632	#endif
2633	iterator1 end () const {
2634	self_type &m = (*this) ();
2635	return m.find1 (1, m.size1 (), index2 ());
2636	}
2637	BOOST_UBLAS_INLINE
2638	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2639	typename self_type::
2640	#endif
2641	reverse_iterator1 rbegin () const {
2642	return reverse_iterator1 (end ());
2643	}
2644	BOOST_UBLAS_INLINE
2645	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2646	typename self_type::
2647	#endif
2648	reverse_iterator1 rend () const {
2649	return reverse_iterator1 (begin ());
2650	}
2651	#endif
2652
2653	// Indices
2654	BOOST_UBLAS_INLINE
2655	size_type index1 () const {
2656	if (rank_ == 1) {
2657	BOOST_UBLAS_CHECK (layout_type::index_M ((*itv_).first, (*it_).first) < (*this) ().size1 (), bad_index ());
2658	return layout_type::index_M ((*itv_).first, (*it_).first);
2659	} else {
2660	return i_;
2661	}
2662	}
2663	BOOST_UBLAS_INLINE
2664	size_type index2 () const {
2665	BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
2666	if (rank_ == 1) {
2667	BOOST_UBLAS_CHECK (layout_type::index_m ((*itv_).first, (*it_).first) < (*this) ().size2 (), bad_index ());
2668	return layout_type::index_m ((*itv_).first, (*it_).first);
2669	} else {
2670	return j_;
2671	}
2672	}
2673
2674	// Assignment
2675	BOOST_UBLAS_INLINE
2676	iterator2 &operator = (const iterator2 &it) {
2677	container_reference<self_type>::assign (&it ());
2678	rank_ = it.rank_;
2679	i_ = it.i_;
2680	j_ = it.j_;
2681	itv_ = it.itv_;
2682	it_ = it.it_;
2683	return *this;
2684	}
2685
2686	// Comparison
2687	BOOST_UBLAS_INLINE
2688	bool operator == (const iterator2 &it) const {
2689	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2690	// BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
2691	if (rank_ == 1 || it.rank_ == 1) {
2692	return it_ == it.it_;
2693	} else {
2694	return i_ == it.i_ && j_ == it.j_;
2695	}
2696	}
2697
2698	private:
2699	int rank_;
2700	size_type i_;
2701	size_type j_;
2702	vector_subiterator_type itv_;
2703	subiterator_type it_;
2704
2705	friend class const_iterator2;
2706	};
2707
2708	BOOST_UBLAS_INLINE
2709	iterator2 begin2 () {
2710	return find2 (0, 0, 0);
2711	}
2712	BOOST_UBLAS_INLINE
2713	iterator2 end2 () {
2714	return find2 (0, 0, size2_);
2715	}
2716
2717	// Reverse iterators
2718
2719	BOOST_UBLAS_INLINE
2720	const_reverse_iterator1 rbegin1 () const {
2721	return const_reverse_iterator1 (end1 ());
2722	}
2723	BOOST_UBLAS_INLINE
2724	const_reverse_iterator1 crbegin1 () const {
2725	return rbegin1 ();
2726	}
2727	BOOST_UBLAS_INLINE
2728	const_reverse_iterator1 rend1 () const {
2729	return const_reverse_iterator1 (begin1 ());
2730	}
2731	BOOST_UBLAS_INLINE
2732	const_reverse_iterator1 crend1 () const {
2733	return rend1 ();
2734	}
2735
2736	BOOST_UBLAS_INLINE
2737	reverse_iterator1 rbegin1 () {
2738	return reverse_iterator1 (end1 ());
2739	}
2740	BOOST_UBLAS_INLINE
2741	reverse_iterator1 rend1 () {
2742	return reverse_iterator1 (begin1 ());
2743	}
2744
2745	BOOST_UBLAS_INLINE
2746	const_reverse_iterator2 rbegin2 () const {
2747	return const_reverse_iterator2 (end2 ());
2748	}
2749	BOOST_UBLAS_INLINE
2750	const_reverse_iterator2 crbegin2 () const {
2751	return rbegin2 ();
2752	}
2753	BOOST_UBLAS_INLINE
2754	const_reverse_iterator2 rend2 () const {
2755	return const_reverse_iterator2 (begin2 ());
2756	}
2757	BOOST_UBLAS_INLINE
2758	const_reverse_iterator2 crend2 () const {
2759	return rend2 ();
2760	}
2761
2762	BOOST_UBLAS_INLINE
2763	reverse_iterator2 rbegin2 () {
2764	return reverse_iterator2 (end2 ());
2765	}
2766	BOOST_UBLAS_INLINE
2767	reverse_iterator2 rend2 () {
2768	return reverse_iterator2 (begin2 ());
2769	}
2770
2771	// Serialization
2772	template<class Archive>
2773	void serialize(Archive & ar, const unsigned int /* file_version */){
2774	serialization::collection_size_type s1 (size1_);
2775	serialization::collection_size_type s2 (size2_);
2776	ar & serialization::make_nvp(n: "size1",v&: s1);
2777	ar & serialization::make_nvp(n: "size2",v&: s2);
2778	if (Archive::is_loading::value) {
2779	size1_ = s1;
2780	size2_ = s2;
2781	}
2782	ar & serialization::make_nvp("data", data_);
2783	}
2784
2785	private:
2786	size_type size1_;
2787	size_type size2_;
2788	array_type data_;
2789	static const value_type zero_;
2790	};
2791
2792	template<class T, class L, class A>
2793	const typename mapped_vector_of_mapped_vector<T, L, A>::value_type mapped_vector_of_mapped_vector<T, L, A>::zero_ = value_type/*zero*/();
2794
2795
2796	// Comperssed array based sparse matrix class
2797	// Thanks to Kresimir Fresl for extending this to cover different index bases.
2798	template<class T, class L, std::size_t IB, class IA, class TA>
2799	class compressed_matrix:
2800	public matrix_container<compressed_matrix<T, L, IB, IA, TA> > {
2801
2802	typedef T &true_reference;
2803	typedef T *pointer;
2804	typedef const T *const_pointer;
2805	typedef L layout_type;
2806	typedef compressed_matrix<T, L, IB, IA, TA> self_type;
2807	public:
2808	#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
2809	using matrix_container<self_type>::operator ();
2810	#endif
2811	// ISSUE require type consistency check
2812	// is_convertable (IA::size_type, TA::size_type)
2813	typedef typename IA::value_type size_type;
2814	// size_type for the data arrays.
2815	typedef typename IA::size_type array_size_type;
2816	// FIXME difference type for sparse storage iterators should it be in the container?
2817	typedef typename IA::difference_type difference_type;
2818	typedef T value_type;
2819	typedef const T &const_reference;
2820	#ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
2821	typedef T &reference;
2822	#else
2823	typedef sparse_matrix_element<self_type> reference;
2824	#endif
2825	typedef IA index_array_type;
2826	typedef TA value_array_type;
2827	typedef const matrix_reference<const self_type> const_closure_type;
2828	typedef matrix_reference<self_type> closure_type;
2829	typedef compressed_vector<T, IB, IA, TA> vector_temporary_type;
2830	typedef self_type matrix_temporary_type;
2831	typedef sparse_tag storage_category;
2832	typedef typename L::orientation_category orientation_category;
2833
2834	// Construction and destruction
2835	BOOST_UBLAS_INLINE
2836	compressed_matrix ():
2837	matrix_container<self_type> (),
2838	size1_ (0), size2_ (0), capacity_ (restrict_capacity (non_zeros: 0)),
2839	filled1_ (1), filled2_ (0),
2840	index1_data_ (layout_type::size_M (size1_, size2_) + 1), index2_data_ (capacity_), value_data_ (capacity_) {
2841	index1_data_ [filled1_ - 1] = k_based (zero_based_index: filled2_);
2842	storage_invariants ();
2843	}
2844	BOOST_UBLAS_INLINE
2845	compressed_matrix (size_type size1, size_type size2, size_type non_zeros = 0):
2846	matrix_container<self_type> (),
2847	size1_ (size1), size2_ (size2), capacity_ (restrict_capacity (non_zeros)),
2848	filled1_ (1), filled2_ (0),
2849	index1_data_ (layout_type::size_M (size1_, size2_) + 1), index2_data_ (capacity_), value_data_ (capacity_) {
2850	index1_data_ [filled1_ - 1] = k_based (zero_based_index: filled2_);
2851	storage_invariants ();
2852	}
2853	BOOST_UBLAS_INLINE
2854	compressed_matrix (const compressed_matrix &m):
2855	matrix_container<self_type> (),
2856	size1_ (m.size1_), size2_ (m.size2_), capacity_ (m.capacity_),
2857	filled1_ (m.filled1_), filled2_ (m.filled2_),
2858	index1_data_ (m.index1_data_), index2_data_ (m.index2_data_), value_data_ (m.value_data_) {
2859	storage_invariants ();
2860	}
2861
2862	BOOST_UBLAS_INLINE
2863	compressed_matrix (const coordinate_matrix<T, L, IB, IA, TA> &m):
2864	matrix_container<self_type> (),
2865	size1_ (m.size1()), size2_ (m.size2()),
2866	index1_data_ (layout_type::size_M (size1_, size2_) + 1)
2867	{
2868	m.sort();
2869	reserve(non_zeros: m.nnz(), preserve: false);
2870	filled2_ = m.nnz();
2871	const_subiterator_type i_start = m.index1_data().begin();
2872	const_subiterator_type i_end = (i_start + filled2_);
2873	const_subiterator_type i = i_start;
2874	size_type r = 1;
2875	for (; (r < layout_type::size_M (size1_, size2_)) && (i != i_end); ++r) {
2876	i = std::lower_bound(i, i_end, r);
2877	index1_data_[r] = k_based( zero_based_index: i - i_start );
2878	}
2879	filled1_ = r + 1;
2880	std::copy( m.index2_data().begin(), m.index2_data().begin() + filled2_, index2_data_.begin());
2881	std::copy( m.value_data().begin(), m.value_data().begin() + filled2_, value_data_.begin());
2882	index1_data_ [filled1_ - 1] = k_based(zero_based_index: filled2_);
2883	storage_invariants ();
2884	}
2885
2886	template<class AE>
2887	BOOST_UBLAS_INLINE
2888	compressed_matrix (const matrix_expression<AE> &ae, size_type non_zeros = 0):
2889	matrix_container<self_type> (),
2890	size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), capacity_ (restrict_capacity (non_zeros)),
2891	filled1_ (1), filled2_ (0),
2892	index1_data_ (layout_type::size_M (ae ().size1 (), ae ().size2 ()) + 1),
2893	index2_data_ (capacity_), value_data_ (capacity_) {
2894	index1_data_ [filled1_ - 1] = k_based (zero_based_index: filled2_);
2895	storage_invariants ();
2896	matrix_assign<scalar_assign> (*this, ae);
2897	}
2898
2899	// Accessors
2900	BOOST_UBLAS_INLINE
2901	size_type size1 () const {
2902	return size1_;
2903	}
2904	BOOST_UBLAS_INLINE
2905	size_type size2 () const {
2906	return size2_;
2907	}
2908	BOOST_UBLAS_INLINE
2909	size_type nnz_capacity () const {
2910	return capacity_;
2911	}
2912	BOOST_UBLAS_INLINE
2913	size_type nnz () const {
2914	return filled2_;
2915	}
2916
2917	// Storage accessors
2918	BOOST_UBLAS_INLINE
2919	static size_type index_base () {
2920	return IB;
2921	}
2922	BOOST_UBLAS_INLINE
2923	array_size_type filled1 () const {
2924	return filled1_;
2925	}
2926	BOOST_UBLAS_INLINE
2927	array_size_type filled2 () const {
2928	return filled2_;
2929	}
2930	BOOST_UBLAS_INLINE
2931	const index_array_type &index1_data () const {
2932	return index1_data_;
2933	}
2934	BOOST_UBLAS_INLINE
2935	const index_array_type &index2_data () const {
2936	return index2_data_;
2937	}
2938	BOOST_UBLAS_INLINE
2939	const value_array_type &value_data () const {
2940	return value_data_;
2941	}
2942	BOOST_UBLAS_INLINE
2943	void set_filled (const array_size_type& filled1, const array_size_type& filled2) {
2944	filled1_ = filled1;
2945	filled2_ = filled2;
2946	storage_invariants ();
2947	}
2948	BOOST_UBLAS_INLINE
2949	index_array_type &index1_data () {
2950	return index1_data_;
2951	}
2952	BOOST_UBLAS_INLINE
2953	index_array_type &index2_data () {
2954	return index2_data_;
2955	}
2956	BOOST_UBLAS_INLINE
2957	value_array_type &value_data () {
2958	return value_data_;
2959	}
2960	BOOST_UBLAS_INLINE
2961	void complete_index1_data () {
2962	while (filled1_ <= layout_type::size_M (size1_, size2_)) {
2963	this->index1_data_ [filled1_] = k_based (zero_based_index: filled2_);
2964	++ this->filled1_;
2965	}
2966	}
2967
2968	// Resizing
2969	private:
2970	BOOST_UBLAS_INLINE
2971	size_type restrict_capacity (size_type non_zeros) const {
2972	non_zeros = (std::max) (non_zeros, (std::min) (size1_, size2_));
2973	// Guarding against overflow - Thanks to Alexei Novakov for the hint.
2974	// non_zeros = (std::min) (non_zeros, size1_ * size2_);
2975	if (size1_ > 0 && non_zeros / size1_ >= size2_)
2976	non_zeros = size1_ * size2_;
2977	return non_zeros;
2978	}
2979	public:
2980	BOOST_UBLAS_INLINE
2981	void resize (size_type size1, size_type size2, bool preserve = true) {
2982	// FIXME preserve unimplemented
2983	BOOST_UBLAS_CHECK (!preserve, internal_logic ());
2984	size1_ = size1;
2985	size2_ = size2;
2986	capacity_ = restrict_capacity (non_zeros: capacity_);
2987	filled1_ = 1;
2988	filled2_ = 0;
2989	index1_data_.resize (layout_type::size_M (size1_, size2_) + 1);
2990	index2_data_.resize (capacity_);
2991	value_data_.resize (capacity_);
2992	index1_data_ [filled1_ - 1] = k_based (zero_based_index: filled2_);
2993	storage_invariants ();
2994	}
2995
2996	// Reserving
2997	BOOST_UBLAS_INLINE
2998	void reserve (size_type non_zeros, bool preserve = true) {
2999	capacity_ = restrict_capacity (non_zeros);
3000	if (preserve) {
3001	index2_data_.resize (capacity_, size_type ());
3002	value_data_.resize (capacity_, value_type ());
3003	filled2_ = (std::min) (capacity_, filled2_);
3004	}
3005	else {
3006	index2_data_.resize (capacity_);
3007	value_data_.resize (capacity_);
3008	filled1_ = 1;
3009	filled2_ = 0;
3010	index1_data_ [filled1_ - 1] = k_based (zero_based_index: filled2_);
3011	}
3012	storage_invariants ();
3013	}
3014
3015	// Element support
3016	BOOST_UBLAS_INLINE
3017	pointer find_element (size_type i, size_type j) {
3018	return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i, j));
3019	}
3020	BOOST_UBLAS_INLINE
3021	const_pointer find_element (size_type i, size_type j) const {
3022	size_type element1 (layout_type::index_M (i, j));
3023	size_type element2 (layout_type::index_m (i, j));
3024	if (filled1_ <= element1 + 1)
3025	return 0;
3026	vector_const_subiterator_type itv (index1_data_.begin () + element1);
3027	const_subiterator_type it_begin (index2_data_.begin () + zero_based (k_based_index: *itv));
3028	const_subiterator_type it_end (index2_data_.begin () + zero_based (k_based_index: *(itv + 1)));
3029	const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (zero_based_index: element2), std::less<size_type> ()));
3030	if (it == it_end || *it != k_based (zero_based_index: element2))
3031	return 0;
3032	return &value_data_ [it - index2_data_.begin ()];
3033	}
3034
3035	// Element access
3036	BOOST_UBLAS_INLINE
3037	const_reference operator () (size_type i, size_type j) const {
3038	const_pointer p = find_element (i, j);
3039	if (p)
3040	return *p;
3041	else
3042	return zero_;
3043	}
3044	BOOST_UBLAS_INLINE
3045	reference operator () (size_type i, size_type j) {
3046	#ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
3047	size_type element1 (layout_type::index_M (i, j));
3048	size_type element2 (layout_type::index_m (i, j));
3049	if (filled1_ <= element1 + 1)
3050	return insert_element (i, j, value_type/*zero*/());
3051	pointer p = find_element (i, j);
3052	if (p)
3053	return *p;
3054	else
3055	return insert_element (i, j, value_type/*zero*/());
3056	#else
3057	return reference (*this, i, j);
3058	#endif
3059	}
3060
3061	// Element assignment
3062	BOOST_UBLAS_INLINE
3063	true_reference insert_element (size_type i, size_type j, const_reference t) {
3064	BOOST_UBLAS_CHECK (!find_element (i, j), bad_index ()); // duplicate element
3065	if (filled2_ >= capacity_)
3066	reserve (non_zeros: 2 * filled2_, preserve: true);
3067	BOOST_UBLAS_CHECK (filled2_ < capacity_, internal_logic ());
3068	size_type element1 = layout_type::index_M (i, j);
3069	size_type element2 = layout_type::index_m (i, j);
3070	while (filled1_ <= element1 + 1) {
3071	index1_data_ [filled1_] = k_based (zero_based_index: filled2_);
3072	++ filled1_;
3073	}
3074	vector_subiterator_type itv (index1_data_.begin () + element1);
3075	subiterator_type it_begin (index2_data_.begin () + zero_based (k_based_index: *itv));
3076	subiterator_type it_end (index2_data_.begin () + zero_based (k_based_index: *(itv + 1)));
3077	subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (zero_based_index: element2), std::less<size_type> ()));
3078	typename std::iterator_traits<subiterator_type>::difference_type n = it - index2_data_.begin ();
3079	BOOST_UBLAS_CHECK (it == it_end || *it != k_based (element2), internal_logic ()); // duplicate bound by lower_bound
3080	++ filled2_;
3081	it = index2_data_.begin () + n;
3082	std::copy_backward (it, index2_data_.begin () + filled2_ - 1, index2_data_.begin () + filled2_);
3083	*it = k_based (zero_based_index: element2);
3084	typename value_array_type::iterator itt (value_data_.begin () + n);
3085	std::copy_backward (itt, value_data_.begin () + filled2_ - 1, value_data_.begin () + filled2_);
3086	*itt = t;
3087	while (element1 + 1 < filled1_) {
3088	++ index1_data_ [element1 + 1];
3089	++ element1;
3090	}
3091	storage_invariants ();
3092	return *itt;
3093	}
3094	BOOST_UBLAS_INLINE
3095	void erase_element (size_type i, size_type j) {
3096	size_type element1 = layout_type::index_M (i, j);
3097	size_type element2 = layout_type::index_m (i, j);
3098	if (element1 + 1 >= filled1_)
3099	return;
3100	vector_subiterator_type itv (index1_data_.begin () + element1);
3101	subiterator_type it_begin (index2_data_.begin () + zero_based (k_based_index: *itv));
3102	subiterator_type it_end (index2_data_.begin () + zero_based (k_based_index: *(itv + 1)));
3103	subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (zero_based_index: element2), std::less<size_type> ()));
3104	if (it != it_end && *it == k_based (zero_based_index: element2)) {
3105	typename std::iterator_traits<subiterator_type>::difference_type n = it - index2_data_.begin ();
3106	std::copy (it + 1, index2_data_.begin () + filled2_, it);
3107	typename value_array_type::iterator itt (value_data_.begin () + n);
3108	std::copy (itt + 1, value_data_.begin () + filled2_, itt);
3109	-- filled2_;
3110	while (index1_data_ [filled1_ - 2] > k_based (zero_based_index: filled2_)) {
3111	index1_data_ [filled1_ - 1] = 0;
3112	-- filled1_;
3113	}
3114	while (element1 + 1 < filled1_) {
3115	-- index1_data_ [element1 + 1];
3116	++ element1;
3117	}
3118	}
3119	storage_invariants ();
3120	}
3121
3122	// Zeroing
3123	BOOST_UBLAS_INLINE
3124	void clear () {
3125	filled1_ = 1;
3126	filled2_ = 0;
3127	index1_data_ [filled1_ - 1] = k_based (zero_based_index: filled2_);
3128	storage_invariants ();
3129	}
3130
3131	// Assignment
3132	BOOST_UBLAS_INLINE
3133	compressed_matrix &operator = (const compressed_matrix &m) {
3134	if (this != &m) {
3135	size1_ = m.size1_;
3136	size2_ = m.size2_;
3137	capacity_ = m.capacity_;
3138	filled1_ = m.filled1_;
3139	filled2_ = m.filled2_;
3140	index1_data_ = m.index1_data_;
3141	index2_data_ = m.index2_data_;
3142	value_data_ = m.value_data_;
3143	}
3144	storage_invariants ();
3145	return *this;
3146	}
3147	template<class C> // Container assignment without temporary
3148	BOOST_UBLAS_INLINE
3149	compressed_matrix &operator = (const matrix_container<C> &m) {
3150	resize (size1: m ().size1 (), size2: m ().size2 (), preserve: false);
3151	assign (m);
3152	return *this;
3153	}
3154	BOOST_UBLAS_INLINE
3155	compressed_matrix &assign_temporary (compressed_matrix &m) {
3156	swap (m);
3157	return *this;
3158	}
3159	template<class AE>
3160	BOOST_UBLAS_INLINE
3161	compressed_matrix &operator = (const matrix_expression<AE> &ae) {
3162	self_type temporary (ae, capacity_);
3163	return assign_temporary (m&: temporary);
3164	}
3165	template<class AE>
3166	BOOST_UBLAS_INLINE
3167	compressed_matrix &assign (const matrix_expression<AE> &ae) {
3168	matrix_assign<scalar_assign> (*this, ae);
3169	return *this;
3170	}
3171	template<class AE>
3172	BOOST_UBLAS_INLINE
3173	compressed_matrix& operator += (const matrix_expression<AE> &ae) {
3174	self_type temporary (*this + ae, capacity_);
3175	return assign_temporary (m&: temporary);
3176	}
3177	template<class C> // Container assignment without temporary
3178	BOOST_UBLAS_INLINE
3179	compressed_matrix &operator += (const matrix_container<C> &m) {
3180	plus_assign (m);
3181	return *this;
3182	}
3183	template<class AE>
3184	BOOST_UBLAS_INLINE
3185	compressed_matrix &plus_assign (const matrix_expression<AE> &ae) {
3186	matrix_assign<scalar_plus_assign> (*this, ae);
3187	return *this;
3188	}
3189	template<class AE>
3190	BOOST_UBLAS_INLINE
3191	compressed_matrix& operator -= (const matrix_expression<AE> &ae) {
3192	self_type temporary (*this - ae, capacity_);
3193	return assign_temporary (m&: temporary);
3194	}
3195	template<class C> // Container assignment without temporary
3196	BOOST_UBLAS_INLINE
3197	compressed_matrix &operator -= (const matrix_container<C> &m) {
3198	minus_assign (m);
3199	return *this;
3200	}
3201	template<class AE>
3202	BOOST_UBLAS_INLINE
3203	compressed_matrix &minus_assign (const matrix_expression<AE> &ae) {
3204	matrix_assign<scalar_minus_assign> (*this, ae);
3205	return *this;
3206	}
3207	template<class AT>
3208	BOOST_UBLAS_INLINE
3209	compressed_matrix& operator *= (const AT &at) {
3210	matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
3211	return *this;
3212	}
3213	template<class AT>
3214	BOOST_UBLAS_INLINE
3215	compressed_matrix& operator /= (const AT &at) {
3216	matrix_assign_scalar<scalar_divides_assign> (*this, at);
3217	return *this;
3218	}
3219
3220	// Swapping
3221	BOOST_UBLAS_INLINE
3222	void swap (compressed_matrix &m) {
3223	if (this != &m) {
3224	std::swap (size1_, m.size1_);
3225	std::swap (size2_, m.size2_);
3226	std::swap (capacity_, m.capacity_);
3227	std::swap (filled1_, m.filled1_);
3228	std::swap (filled2_, m.filled2_);
3229	index1_data_.swap (m.index1_data_);
3230	index2_data_.swap (m.index2_data_);
3231	value_data_.swap (m.value_data_);
3232	}
3233	storage_invariants ();
3234	}
3235	BOOST_UBLAS_INLINE
3236	friend void swap (compressed_matrix &m1, compressed_matrix &m2) {
3237	m1.swap (m2);
3238	}
3239
3240	// Back element insertion and erasure
3241	BOOST_UBLAS_INLINE
3242	void push_back (size_type i, size_type j, const_reference t) {
3243	if (filled2_ >= capacity_)
3244	reserve (non_zeros: 2 * filled2_, preserve: true);
3245	BOOST_UBLAS_CHECK (filled2_ < capacity_, internal_logic ());
3246	size_type element1 = layout_type::index_M (i, j);
3247	size_type element2 = layout_type::index_m (i, j);
3248	while (filled1_ < element1 + 2) {
3249	index1_data_ [filled1_] = k_based (zero_based_index: filled2_);
3250	++ filled1_;
3251	}
3252	// must maintain sort order
3253	BOOST_UBLAS_CHECK ((filled1_ == element1 + 2 &&
3254	(filled2_ == zero_based (index1_data_ [filled1_ - 2]) ||
3255	index2_data_ [filled2_ - 1] < k_based (element2))), external_logic ());
3256	++ filled2_;
3257	index1_data_ [filled1_ - 1] = k_based (zero_based_index: filled2_);
3258	index2_data_ [filled2_ - 1] = k_based (zero_based_index: element2);
3259	value_data_ [filled2_ - 1] = t;
3260	storage_invariants ();
3261	}
3262	BOOST_UBLAS_INLINE
3263	void pop_back () {
3264	BOOST_UBLAS_CHECK (filled1_ > 0 && filled2_ > 0, external_logic ());
3265	-- filled2_;
3266	while (index1_data_ [filled1_ - 2] > k_based (zero_based_index: filled2_)) {
3267	index1_data_ [filled1_ - 1] = 0;
3268	-- filled1_;
3269	}
3270	-- index1_data_ [filled1_ - 1];
3271	storage_invariants ();
3272	}
3273
3274	// Iterator types
3275	private:
3276	// Use index array iterator
3277	typedef typename IA::const_iterator vector_const_subiterator_type;
3278	typedef typename IA::iterator vector_subiterator_type;
3279	typedef typename IA::const_iterator const_subiterator_type;
3280	typedef typename IA::iterator subiterator_type;
3281
3282	BOOST_UBLAS_INLINE
3283	true_reference at_element (size_type i, size_type j) {
3284	pointer p = find_element (i, j);
3285	BOOST_UBLAS_CHECK (p, bad_index ());
3286	return *p;
3287	}
3288
3289	public:
3290	class const_iterator1;
3291	class iterator1;
3292	class const_iterator2;
3293	class iterator2;
3294	typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
3295	typedef reverse_iterator_base1<iterator1> reverse_iterator1;
3296	typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
3297	typedef reverse_iterator_base2<iterator2> reverse_iterator2;
3298
3299	// Element lookup
3300	// BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
3301	const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const {
3302	for (;;) {
3303	array_size_type address1 (layout_type::index_M (i, j));
3304	array_size_type address2 (layout_type::index_m (i, j));
3305	vector_const_subiterator_type itv (index1_data_.begin () + (std::min) (filled1_ - 1, address1));
3306	if (filled1_ <= address1 + 1)
3307	return const_iterator1 (*this, rank, i, j, itv, index2_data_.begin () + filled2_);
3308
3309	const_subiterator_type it_begin (index2_data_.begin () + zero_based (k_based_index: *itv));
3310	const_subiterator_type it_end (index2_data_.begin () + zero_based (k_based_index: *(itv + 1)));
3311
3312	const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (zero_based_index: address2), std::less<size_type> ()));
3313	if (rank == 0)
3314	return const_iterator1 (*this, rank, i, j, itv, it);
3315	if (it != it_end && zero_based (k_based_index: *it) == address2)
3316	return const_iterator1 (*this, rank, i, j, itv, it);
3317	if (direction > 0) {
3318	if (layout_type::fast_i ()) {
3319	if (it == it_end)
3320	return const_iterator1 (*this, rank, i, j, itv, it);
3321	i = zero_based (k_based_index: *it);
3322	} else {
3323	if (i >= size1_)
3324	return const_iterator1 (*this, rank, i, j, itv, it);
3325	++ i;
3326	}
3327	} else /* if (direction < 0) */ {
3328	if (layout_type::fast_i ()) {
3329	if (it == index2_data_.begin () + zero_based (k_based_index: *itv))
3330	return const_iterator1 (*this, rank, i, j, itv, it);
3331	i = zero_based (k_based_index: *(it - 1));
3332	} else {
3333	if (i == 0)
3334	return const_iterator1 (*this, rank, i, j, itv, it);
3335	-- i;
3336	}
3337	}
3338	}
3339	}
3340	// BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
3341	iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) {
3342	for (;;) {
3343	array_size_type address1 (layout_type::index_M (i, j));
3344	array_size_type address2 (layout_type::index_m (i, j));
3345	vector_subiterator_type itv (index1_data_.begin () + (std::min) (filled1_ - 1, address1));
3346	if (filled1_ <= address1 + 1)
3347	return iterator1 (*this, rank, i, j, itv, index2_data_.begin () + filled2_);
3348
3349	subiterator_type it_begin (index2_data_.begin () + zero_based (k_based_index: *itv));
3350	subiterator_type it_end (index2_data_.begin () + zero_based (k_based_index: *(itv + 1)));
3351
3352	subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (zero_based_index: address2), std::less<size_type> ()));
3353	if (rank == 0)
3354	return iterator1 (*this, rank, i, j, itv, it);
3355	if (it != it_end && zero_based (k_based_index: *it) == address2)
3356	return iterator1 (*this, rank, i, j, itv, it);
3357	if (direction > 0) {
3358	if (layout_type::fast_i ()) {
3359	if (it == it_end)
3360	return iterator1 (*this, rank, i, j, itv, it);
3361	i = zero_based (k_based_index: *it);
3362	} else {
3363	if (i >= size1_)
3364	return iterator1 (*this, rank, i, j, itv, it);
3365	++ i;
3366	}
3367	} else /* if (direction < 0) */ {
3368	if (layout_type::fast_i ()) {
3369	if (it == index2_data_.begin () + zero_based (k_based_index: *itv))
3370	return iterator1 (*this, rank, i, j, itv, it);
3371	i = zero_based (k_based_index: *(it - 1));
3372	} else {
3373	if (i == 0)
3374	return iterator1 (*this, rank, i, j, itv, it);
3375	-- i;
3376	}
3377	}
3378	}
3379	}
3380	// BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
3381	const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const {
3382	for (;;) {
3383	array_size_type address1 (layout_type::index_M (i, j));
3384	array_size_type address2 (layout_type::index_m (i, j));
3385	vector_const_subiterator_type itv (index1_data_.begin () + (std::min) (filled1_ - 1, address1));
3386	if (filled1_ <= address1 + 1)
3387	return const_iterator2 (*this, rank, i, j, itv, index2_data_.begin () + filled2_);
3388
3389	const_subiterator_type it_begin (index2_data_.begin () + zero_based (k_based_index: *itv));
3390	const_subiterator_type it_end (index2_data_.begin () + zero_based (k_based_index: *(itv + 1)));
3391
3392	const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (zero_based_index: address2), std::less<size_type> ()));
3393	if (rank == 0)
3394	return const_iterator2 (*this, rank, i, j, itv, it);
3395	if (it != it_end && zero_based (k_based_index: *it) == address2)
3396	return const_iterator2 (*this, rank, i, j, itv, it);
3397	if (direction > 0) {
3398	if (layout_type::fast_j ()) {
3399	if (it == it_end)
3400	return const_iterator2 (*this, rank, i, j, itv, it);
3401	j = zero_based (k_based_index: *it);
3402	} else {
3403	if (j >= size2_)
3404	return const_iterator2 (*this, rank, i, j, itv, it);
3405	++ j;
3406	}
3407	} else /* if (direction < 0) */ {
3408	if (layout_type::fast_j ()) {
3409	if (it == index2_data_.begin () + zero_based (k_based_index: *itv))
3410	return const_iterator2 (*this, rank, i, j, itv, it);
3411	j = zero_based (k_based_index: *(it - 1));
3412	} else {
3413	if (j == 0)
3414	return const_iterator2 (*this, rank, i, j, itv, it);
3415	-- j;
3416	}
3417	}
3418	}
3419	}
3420	// BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
3421	iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) {
3422	for (;;) {
3423	array_size_type address1 (layout_type::index_M (i, j));
3424	array_size_type address2 (layout_type::index_m (i, j));
3425	vector_subiterator_type itv (index1_data_.begin () + (std::min) (filled1_ - 1, address1));
3426	if (filled1_ <= address1 + 1)
3427	return iterator2 (*this, rank, i, j, itv, index2_data_.begin () + filled2_);
3428
3429	subiterator_type it_begin (index2_data_.begin () + zero_based (k_based_index: *itv));
3430	subiterator_type it_end (index2_data_.begin () + zero_based (k_based_index: *(itv + 1)));
3431
3432	subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (zero_based_index: address2), std::less<size_type> ()));
3433	if (rank == 0)
3434	return iterator2 (*this, rank, i, j, itv, it);
3435	if (it != it_end && zero_based (k_based_index: *it) == address2)
3436	return iterator2 (*this, rank, i, j, itv, it);
3437	if (direction > 0) {
3438	if (layout_type::fast_j ()) {
3439	if (it == it_end)
3440	return iterator2 (*this, rank, i, j, itv, it);
3441	j = zero_based (k_based_index: *it);
3442	} else {
3443	if (j >= size2_)
3444	return iterator2 (*this, rank, i, j, itv, it);
3445	++ j;
3446	}
3447	} else /* if (direction < 0) */ {
3448	if (layout_type::fast_j ()) {
3449	if (it == index2_data_.begin () + zero_based (k_based_index: *itv))
3450	return iterator2 (*this, rank, i, j, itv, it);
3451	j = zero_based (k_based_index: *(it - 1));
3452	} else {
3453	if (j == 0)
3454	return iterator2 (*this, rank, i, j, itv, it);
3455	-- j;
3456	}
3457	}
3458	}
3459	}
3460
3461
3462	class const_iterator1:
3463	public container_const_reference<compressed_matrix>,
3464	public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
3465	const_iterator1, value_type> {
3466	public:
3467	typedef typename compressed_matrix::value_type value_type;
3468	typedef typename compressed_matrix::difference_type difference_type;
3469	typedef typename compressed_matrix::const_reference reference;
3470	typedef const typename compressed_matrix::pointer pointer;
3471
3472	typedef const_iterator2 dual_iterator_type;
3473	typedef const_reverse_iterator2 dual_reverse_iterator_type;
3474
3475	// Construction and destruction
3476	BOOST_UBLAS_INLINE
3477	const_iterator1 ():
3478	container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
3479	BOOST_UBLAS_INLINE
3480	const_iterator1 (const self_type &m, int rank, size_type i, size_type j, const vector_const_subiterator_type &itv, const const_subiterator_type &it):
3481	container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
3482	BOOST_UBLAS_INLINE
3483	const_iterator1 (const iterator1 &it):
3484	container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {}
3485
3486	// Arithmetic
3487	BOOST_UBLAS_INLINE
3488	const_iterator1 &operator ++ () {
3489	if (rank_ == 1 && layout_type::fast_i ())
3490	++ it_;
3491	else {
3492	i_ = index1 () + 1;
3493	if (rank_ == 1)
3494	*this = (*this) ().find1 (rank_, i_, j_, 1);
3495	}
3496	return *this;
3497	}
3498	BOOST_UBLAS_INLINE
3499	const_iterator1 &operator -- () {
3500	if (rank_ == 1 && layout_type::fast_i ())
3501	-- it_;
3502	else {
3503	--i_;
3504	if (rank_ == 1)
3505	*this = (*this) ().find1 (rank_, i_, j_, -1);
3506	}
3507	return *this;
3508	}
3509
3510	// Dereference
3511	BOOST_UBLAS_INLINE
3512	const_reference operator * () const {
3513	BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
3514	BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
3515	if (rank_ == 1) {
3516	return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
3517	} else {
3518	return (*this) () (i_, j_);
3519	}
3520	}
3521
3522	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
3523	BOOST_UBLAS_INLINE
3524	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3525	typename self_type::
3526	#endif
3527	const_iterator2 begin () const {
3528	const self_type &m = (*this) ();
3529	return m.find2 (1, index1 (), 0);
3530	}
3531	BOOST_UBLAS_INLINE
3532	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3533	typename self_type::
3534	#endif
3535	const_iterator2 cbegin () const {
3536	return begin ();
3537	}
3538	BOOST_UBLAS_INLINE
3539	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3540	typename self_type::
3541	#endif
3542	const_iterator2 end () const {
3543	const self_type &m = (*this) ();
3544	return m.find2 (1, index1 (), m.size2 ());
3545	}
3546	BOOST_UBLAS_INLINE
3547	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3548	typename self_type::
3549	#endif
3550	const_iterator2 cend () const {
3551	return end ();
3552	}
3553	BOOST_UBLAS_INLINE
3554	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3555	typename self_type::
3556	#endif
3557	const_reverse_iterator2 rbegin () const {
3558	return const_reverse_iterator2 (end ());
3559	}
3560	BOOST_UBLAS_INLINE
3561	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3562	typename self_type::
3563	#endif
3564	const_reverse_iterator2 crbegin () const {
3565	return rbegin ();
3566	}
3567	BOOST_UBLAS_INLINE
3568	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3569	typename self_type::
3570	#endif
3571	const_reverse_iterator2 rend () const {
3572	return const_reverse_iterator2 (begin ());
3573	}
3574	BOOST_UBLAS_INLINE
3575	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3576	typename self_type::
3577	#endif
3578	const_reverse_iterator2 crend () const {
3579	return rend ();
3580	}
3581	#endif
3582
3583	// Indices
3584	BOOST_UBLAS_INLINE
3585	size_type index1 () const {
3586	BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
3587	if (rank_ == 1) {
3588	BOOST_UBLAS_CHECK (layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
3589	return layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
3590	} else {
3591	return i_;
3592	}
3593	}
3594	BOOST_UBLAS_INLINE
3595	size_type index2 () const {
3596	if (rank_ == 1) {
3597	BOOST_UBLAS_CHECK (layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
3598	return layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
3599	} else {
3600	return j_;
3601	}
3602	}
3603
3604	// Assignment
3605	BOOST_UBLAS_INLINE
3606	const_iterator1 &operator = (const const_iterator1 &it) {
3607	container_const_reference<self_type>::assign (&it ());
3608	rank_ = it.rank_;
3609	i_ = it.i_;
3610	j_ = it.j_;
3611	itv_ = it.itv_;
3612	it_ = it.it_;
3613	return *this;
3614	}
3615
3616	// Comparison
3617	BOOST_UBLAS_INLINE
3618	bool operator == (const const_iterator1 &it) const {
3619	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
3620	// BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
3621	if (rank_ == 1 || it.rank_ == 1) {
3622	return it_ == it.it_;
3623	} else {
3624	return i_ == it.i_ && j_ == it.j_;
3625	}
3626	}
3627
3628	private:
3629	int rank_;
3630	size_type i_;
3631	size_type j_;
3632	vector_const_subiterator_type itv_;
3633	const_subiterator_type it_;
3634	};
3635
3636	BOOST_UBLAS_INLINE
3637	const_iterator1 begin1 () const {
3638	return find1 (0, 0, 0);
3639	}
3640	BOOST_UBLAS_INLINE
3641	const_iterator1 cbegin1 () const {
3642	return begin1 ();
3643	}
3644	BOOST_UBLAS_INLINE
3645	const_iterator1 end1 () const {
3646	return find1 (0, size1_, 0);
3647	}
3648	BOOST_UBLAS_INLINE
3649	const_iterator1 cend1 () const {
3650	return end1 ();
3651	}
3652
3653	class iterator1:
3654	public container_reference<compressed_matrix>,
3655	public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
3656	iterator1, value_type> {
3657	public:
3658	typedef typename compressed_matrix::value_type value_type;
3659	typedef typename compressed_matrix::difference_type difference_type;
3660	typedef typename compressed_matrix::true_reference reference;
3661	typedef typename compressed_matrix::pointer pointer;
3662
3663	typedef iterator2 dual_iterator_type;
3664	typedef reverse_iterator2 dual_reverse_iterator_type;
3665
3666	// Construction and destruction
3667	BOOST_UBLAS_INLINE
3668	iterator1 ():
3669	container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
3670	BOOST_UBLAS_INLINE
3671	iterator1 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it):
3672	container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
3673
3674	// Arithmetic
3675	BOOST_UBLAS_INLINE
3676	iterator1 &operator ++ () {
3677	if (rank_ == 1 && layout_type::fast_i ())
3678	++ it_;
3679	else {
3680	i_ = index1 () + 1;
3681	if (rank_ == 1)
3682	*this = (*this) ().find1 (rank_, i_, j_, 1);
3683	}
3684	return *this;
3685	}
3686	BOOST_UBLAS_INLINE
3687	iterator1 &operator -- () {
3688	if (rank_ == 1 && layout_type::fast_i ())
3689	-- it_;
3690	else {
3691	--i_;
3692	if (rank_ == 1)
3693	*this = (*this) ().find1 (rank_, i_, j_, -1);
3694	}
3695	return *this;
3696	}
3697
3698	// Dereference
3699	BOOST_UBLAS_INLINE
3700	reference operator * () const {
3701	BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
3702	BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
3703	if (rank_ == 1) {
3704	return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
3705	} else {
3706	return (*this) ().at_element (i_, j_);
3707	}
3708	}
3709
3710	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
3711	BOOST_UBLAS_INLINE
3712	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3713	typename self_type::
3714	#endif
3715	iterator2 begin () const {
3716	self_type &m = (*this) ();
3717	return m.find2 (1, index1 (), 0);
3718	}
3719	BOOST_UBLAS_INLINE
3720	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3721	typename self_type::
3722	#endif
3723	iterator2 end () const {
3724	self_type &m = (*this) ();
3725	return m.find2 (1, index1 (), m.size2 ());
3726	}
3727	BOOST_UBLAS_INLINE
3728	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3729	typename self_type::
3730	#endif
3731	reverse_iterator2 rbegin () const {
3732	return reverse_iterator2 (end ());
3733	}
3734	BOOST_UBLAS_INLINE
3735	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3736	typename self_type::
3737	#endif
3738	reverse_iterator2 rend () const {
3739	return reverse_iterator2 (begin ());
3740	}
3741	#endif
3742
3743	// Indices
3744	BOOST_UBLAS_INLINE
3745	size_type index1 () const {
3746	BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
3747	if (rank_ == 1) {
3748	BOOST_UBLAS_CHECK (layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
3749	return layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
3750	} else {
3751	return i_;
3752	}
3753	}
3754	BOOST_UBLAS_INLINE
3755	size_type index2 () const {
3756	if (rank_ == 1) {
3757	BOOST_UBLAS_CHECK (layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
3758	return layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
3759	} else {
3760	return j_;
3761	}
3762	}
3763
3764	// Assignment
3765	BOOST_UBLAS_INLINE
3766	iterator1 &operator = (const iterator1 &it) {
3767	container_reference<self_type>::assign (&it ());
3768	rank_ = it.rank_;
3769	i_ = it.i_;
3770	j_ = it.j_;
3771	itv_ = it.itv_;
3772	it_ = it.it_;
3773	return *this;
3774	}
3775
3776	// Comparison
3777	BOOST_UBLAS_INLINE
3778	bool operator == (const iterator1 &it) const {
3779	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
3780	// BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
3781	if (rank_ == 1 || it.rank_ == 1) {
3782	return it_ == it.it_;
3783	} else {
3784	return i_ == it.i_ && j_ == it.j_;
3785	}
3786	}
3787
3788	private:
3789	int rank_;
3790	size_type i_;
3791	size_type j_;
3792	vector_subiterator_type itv_;
3793	subiterator_type it_;
3794
3795	friend class const_iterator1;
3796	};
3797
3798	BOOST_UBLAS_INLINE
3799	iterator1 begin1 () {
3800	return find1 (0, 0, 0);
3801	}
3802	BOOST_UBLAS_INLINE
3803	iterator1 end1 () {
3804	return find1 (0, size1_, 0);
3805	}
3806
3807	class const_iterator2:
3808	public container_const_reference<compressed_matrix>,
3809	public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
3810	const_iterator2, value_type> {
3811	public:
3812	typedef typename compressed_matrix::value_type value_type;
3813	typedef typename compressed_matrix::difference_type difference_type;
3814	typedef typename compressed_matrix::const_reference reference;
3815	typedef const typename compressed_matrix::pointer pointer;
3816
3817	typedef const_iterator1 dual_iterator_type;
3818	typedef const_reverse_iterator1 dual_reverse_iterator_type;
3819
3820	// Construction and destruction
3821	BOOST_UBLAS_INLINE
3822	const_iterator2 ():
3823	container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
3824	BOOST_UBLAS_INLINE
3825	const_iterator2 (const self_type &m, int rank, size_type i, size_type j, const vector_const_subiterator_type itv, const const_subiterator_type &it):
3826	container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
3827	BOOST_UBLAS_INLINE
3828	const_iterator2 (const iterator2 &it):
3829	container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {}
3830
3831	// Arithmetic
3832	BOOST_UBLAS_INLINE
3833	const_iterator2 &operator ++ () {
3834	if (rank_ == 1 && layout_type::fast_j ())
3835	++ it_;
3836	else {
3837	j_ = index2 () + 1;
3838	if (rank_ == 1)
3839	*this = (*this) ().find2 (rank_, i_, j_, 1);
3840	}
3841	return *this;
3842	}
3843	BOOST_UBLAS_INLINE
3844	const_iterator2 &operator -- () {
3845	if (rank_ == 1 && layout_type::fast_j ())
3846	-- it_;
3847	else {
3848	--j_;
3849	if (rank_ == 1)
3850	*this = (*this) ().find2 (rank_, i_, j_, -1);
3851	}
3852	return *this;
3853	}
3854
3855	// Dereference
3856	BOOST_UBLAS_INLINE
3857	const_reference operator * () const {
3858	BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
3859	BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
3860	if (rank_ == 1) {
3861	return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
3862	} else {
3863	return (*this) () (i_, j_);
3864	}
3865	}
3866
3867	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
3868	BOOST_UBLAS_INLINE
3869	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3870	typename self_type::
3871	#endif
3872	const_iterator1 begin () const {
3873	const self_type &m = (*this) ();
3874	return m.find1 (1, 0, index2 ());
3875	}
3876	BOOST_UBLAS_INLINE
3877	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3878	typename self_type::
3879	#endif
3880	const_iterator1 cbegin () const {
3881	return begin ();
3882	}
3883	BOOST_UBLAS_INLINE
3884	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3885	typename self_type::
3886	#endif
3887	const_iterator1 end () const {
3888	const self_type &m = (*this) ();
3889	return m.find1 (1, m.size1 (), index2 ());
3890	}
3891	BOOST_UBLAS_INLINE
3892	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3893	typename self_type::
3894	#endif
3895	const_iterator1 cend () const {
3896	return end ();
3897	}
3898	BOOST_UBLAS_INLINE
3899	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3900	typename self_type::
3901	#endif
3902	const_reverse_iterator1 rbegin () const {
3903	return const_reverse_iterator1 (end ());
3904	}
3905	BOOST_UBLAS_INLINE
3906	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3907	typename self_type::
3908	#endif
3909	const_reverse_iterator1 crbegin () const {
3910	return rbegin ();
3911	}
3912	BOOST_UBLAS_INLINE
3913	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3914	typename self_type::
3915	#endif
3916	const_reverse_iterator1 rend () const {
3917	return const_reverse_iterator1 (begin ());
3918	}
3919	BOOST_UBLAS_INLINE
3920	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3921	typename self_type::
3922	#endif
3923	const_reverse_iterator1 crend () const {
3924	return rend ();
3925	}
3926	#endif
3927
3928	// Indices
3929	BOOST_UBLAS_INLINE
3930	size_type index1 () const {
3931	if (rank_ == 1) {
3932	BOOST_UBLAS_CHECK (layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
3933	return layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
3934	} else {
3935	return i_;
3936	}
3937	}
3938	BOOST_UBLAS_INLINE
3939	size_type index2 () const {
3940	BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
3941	if (rank_ == 1) {
3942	BOOST_UBLAS_CHECK (layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
3943	return layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
3944	} else {
3945	return j_;
3946	}
3947	}
3948
3949	// Assignment
3950	BOOST_UBLAS_INLINE
3951	const_iterator2 &operator = (const const_iterator2 &it) {
3952	container_const_reference<self_type>::assign (&it ());
3953	rank_ = it.rank_;
3954	i_ = it.i_;
3955	j_ = it.j_;
3956	itv_ = it.itv_;
3957	it_ = it.it_;
3958	return *this;
3959	}
3960
3961	// Comparison
3962	BOOST_UBLAS_INLINE
3963	bool operator == (const const_iterator2 &it) const {
3964	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
3965	// BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
3966	if (rank_ == 1 || it.rank_ == 1) {
3967	return it_ == it.it_;
3968	} else {
3969	return i_ == it.i_ && j_ == it.j_;
3970	}
3971	}
3972
3973	private:
3974	int rank_;
3975	size_type i_;
3976	size_type j_;
3977	vector_const_subiterator_type itv_;
3978	const_subiterator_type it_;
3979	};
3980
3981	BOOST_UBLAS_INLINE
3982	const_iterator2 begin2 () const {
3983	return find2 (0, 0, 0);
3984	}
3985	BOOST_UBLAS_INLINE
3986	const_iterator2 cbegin2 () const {
3987	return begin2 ();
3988	}
3989	BOOST_UBLAS_INLINE
3990	const_iterator2 end2 () const {
3991	return find2 (0, 0, size2_);
3992	}
3993	BOOST_UBLAS_INLINE
3994	const_iterator2 cend2 () const {
3995	return end2 ();
3996	}
3997
3998	class iterator2:
3999	public container_reference<compressed_matrix>,
4000	public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
4001	iterator2, value_type> {
4002	public:
4003	typedef typename compressed_matrix::value_type value_type;
4004	typedef typename compressed_matrix::difference_type difference_type;
4005	typedef typename compressed_matrix::true_reference reference;
4006	typedef typename compressed_matrix::pointer pointer;
4007
4008	typedef iterator1 dual_iterator_type;
4009	typedef reverse_iterator1 dual_reverse_iterator_type;
4010
4011	// Construction and destruction
4012	BOOST_UBLAS_INLINE
4013	iterator2 ():
4014	container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
4015	BOOST_UBLAS_INLINE
4016	iterator2 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it):
4017	container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
4018
4019	// Arithmetic
4020	BOOST_UBLAS_INLINE
4021	iterator2 &operator ++ () {
4022	if (rank_ == 1 && layout_type::fast_j ())
4023	++ it_;
4024	else {
4025	j_ = index2 () + 1;
4026	if (rank_ == 1)
4027	*this = (*this) ().find2 (rank_, i_, j_, 1);
4028	}
4029	return *this;
4030	}
4031	BOOST_UBLAS_INLINE
4032	iterator2 &operator -- () {
4033	if (rank_ == 1 && layout_type::fast_j ())
4034	-- it_;
4035	else {
4036	--j_;
4037	if (rank_ == 1)
4038	*this = (*this) ().find2 (rank_, i_, j_, -1);
4039	}
4040	return *this;
4041	}
4042
4043	// Dereference
4044	BOOST_UBLAS_INLINE
4045	reference operator * () const {
4046	BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
4047	BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
4048	if (rank_ == 1) {
4049	return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
4050	} else {
4051	return (*this) ().at_element (i_, j_);
4052	}
4053	}
4054
4055	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
4056	BOOST_UBLAS_INLINE
4057	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
4058	typename self_type::
4059	#endif
4060	iterator1 begin () const {
4061	self_type &m = (*this) ();
4062	return m.find1 (1, 0, index2 ());
4063	}
4064	BOOST_UBLAS_INLINE
4065	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
4066	typename self_type::
4067	#endif
4068	iterator1 end () const {
4069	self_type &m = (*this) ();
4070	return m.find1 (1, m.size1 (), index2 ());
4071	}
4072	BOOST_UBLAS_INLINE
4073	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
4074	typename self_type::
4075	#endif
4076	reverse_iterator1 rbegin () const {
4077	return reverse_iterator1 (end ());
4078	}
4079	BOOST_UBLAS_INLINE
4080	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
4081	typename self_type::
4082	#endif
4083	reverse_iterator1 rend () const {
4084	return reverse_iterator1 (begin ());
4085	}
4086	#endif
4087
4088	// Indices
4089	BOOST_UBLAS_INLINE
4090	size_type index1 () const {
4091	if (rank_ == 1) {
4092	BOOST_UBLAS_CHECK (layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
4093	return layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
4094	} else {
4095	return i_;
4096	}
4097	}
4098	BOOST_UBLAS_INLINE
4099	size_type index2 () const {
4100	BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
4101	if (rank_ == 1) {
4102	BOOST_UBLAS_CHECK (layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
4103	return layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_));
4104	} else {
4105	return j_;
4106	}
4107	}
4108
4109	// Assignment
4110	BOOST_UBLAS_INLINE
4111	iterator2 &operator = (const iterator2 &it) {
4112	container_reference<self_type>::assign (&it ());
4113	rank_ = it.rank_;
4114	i_ = it.i_;
4115	j_ = it.j_;
4116	itv_ = it.itv_;
4117	it_ = it.it_;
4118	return *this;
4119	}
4120
4121	// Comparison
4122	BOOST_UBLAS_INLINE
4123	bool operator == (const iterator2 &it) const {
4124	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
4125	// BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
4126	if (rank_ == 1 || it.rank_ == 1) {
4127	return it_ == it.it_;
4128	} else {
4129	return i_ == it.i_ && j_ == it.j_;
4130	}
4131	}
4132
4133	private:
4134	int rank_;
4135	size_type i_;
4136	size_type j_;
4137	vector_subiterator_type itv_;
4138	subiterator_type it_;
4139
4140	friend class const_iterator2;
4141	};
4142
4143	BOOST_UBLAS_INLINE
4144	iterator2 begin2 () {
4145	return find2 (0, 0, 0);
4146	}
4147	BOOST_UBLAS_INLINE
4148	iterator2 end2 () {
4149	return find2 (0, 0, size2_);
4150	}
4151
4152	// Reverse iterators
4153
4154	BOOST_UBLAS_INLINE
4155	const_reverse_iterator1 rbegin1 () const {
4156	return const_reverse_iterator1 (end1 ());
4157	}
4158	BOOST_UBLAS_INLINE
4159	const_reverse_iterator1 crbegin1 () const {
4160	return rbegin1 ();
4161	}
4162	BOOST_UBLAS_INLINE
4163	const_reverse_iterator1 rend1 () const {
4164	return const_reverse_iterator1 (begin1 ());
4165	}
4166	BOOST_UBLAS_INLINE
4167	const_reverse_iterator1 crend1 () const {
4168	return rend1 ();
4169	}
4170
4171	BOOST_UBLAS_INLINE
4172	reverse_iterator1 rbegin1 () {
4173	return reverse_iterator1 (end1 ());
4174	}
4175	BOOST_UBLAS_INLINE
4176	reverse_iterator1 rend1 () {
4177	return reverse_iterator1 (begin1 ());
4178	}
4179
4180	BOOST_UBLAS_INLINE
4181	const_reverse_iterator2 rbegin2 () const {
4182	return const_reverse_iterator2 (end2 ());
4183	}
4184	BOOST_UBLAS_INLINE
4185	const_reverse_iterator2 crbegin2 () const {
4186	return rbegin2 ();
4187	}
4188	BOOST_UBLAS_INLINE
4189	const_reverse_iterator2 rend2 () const {
4190	return const_reverse_iterator2 (begin2 ());
4191	}
4192	BOOST_UBLAS_INLINE
4193	const_reverse_iterator2 crend2 () const {
4194	return rend2 ();
4195	}
4196
4197	BOOST_UBLAS_INLINE
4198	reverse_iterator2 rbegin2 () {
4199	return reverse_iterator2 (end2 ());
4200	}
4201	BOOST_UBLAS_INLINE
4202	reverse_iterator2 rend2 () {
4203	return reverse_iterator2 (begin2 ());
4204	}
4205
4206	// Serialization
4207	template<class Archive>
4208	void serialize(Archive & ar, const unsigned int /* file_version */){
4209	serialization::collection_size_type s1 (size1_);
4210	serialization::collection_size_type s2 (size2_);
4211	ar & serialization::make_nvp(n: "size1",v&: s1);
4212	ar & serialization::make_nvp(n: "size2",v&: s2);
4213	if (Archive::is_loading::value) {
4214	size1_ = s1;
4215	size2_ = s2;
4216	}
4217	ar & serialization::make_nvp("capacity", capacity_);
4218	ar & serialization::make_nvp("filled1", filled1_);
4219	ar & serialization::make_nvp("filled2", filled2_);
4220	ar & serialization::make_nvp("index1_data", index1_data_);
4221	ar & serialization::make_nvp("index2_data", index2_data_);
4222	ar & serialization::make_nvp("value_data", value_data_);
4223	storage_invariants();
4224	}
4225
4226	private:
4227	void storage_invariants () const {
4228	BOOST_UBLAS_CHECK (layout_type::size_M (size1_, size2_) + 1 == index1_data_.size (), internal_logic ());
4229	BOOST_UBLAS_CHECK (capacity_ == index2_data_.size (), internal_logic ());
4230	BOOST_UBLAS_CHECK (capacity_ == value_data_.size (), internal_logic ());
4231	BOOST_UBLAS_CHECK (filled1_ > 0 && filled1_ <= layout_type::size_M (size1_, size2_) + 1, internal_logic ());
4232	BOOST_UBLAS_CHECK (filled2_ <= capacity_, internal_logic ());
4233	BOOST_UBLAS_CHECK (index1_data_ [filled1_ - 1] == k_based (filled2_), internal_logic ());
4234	}
4235
4236	size_type size1_;
4237	size_type size2_;
4238	array_size_type capacity_;
4239	array_size_type filled1_;
4240	array_size_type filled2_;
4241	index_array_type index1_data_;
4242	index_array_type index2_data_;
4243	value_array_type value_data_;
4244	static const value_type zero_;
4245
4246	BOOST_UBLAS_INLINE
4247	static size_type zero_based (size_type k_based_index) {
4248	return k_based_index - IB;
4249	}
4250	BOOST_UBLAS_INLINE
4251	static size_type k_based (size_type zero_based_index) {
4252	return zero_based_index + IB;
4253	}
4254
4255	friend class iterator1;
4256	friend class iterator2;
4257	friend class const_iterator1;
4258	friend class const_iterator2;
4259	};
4260
4261	template<class T, class L, std::size_t IB, class IA, class TA>
4262	const typename compressed_matrix<T, L, IB, IA, TA>::value_type compressed_matrix<T, L, IB, IA, TA>::zero_ = value_type/*zero*/();
4263
4264
4265	// Coordinate array based sparse matrix class
4266	// Thanks to Kresimir Fresl for extending this to cover different index bases.
4267	template<class T, class L, std::size_t IB, class IA, class TA>
4268	class coordinate_matrix:
4269	public matrix_container<coordinate_matrix<T, L, IB, IA, TA> > {
4270
4271	typedef T &true_reference;
4272	typedef T *pointer;
4273	typedef const T *const_pointer;
4274	typedef L layout_type;
4275	typedef coordinate_matrix<T, L, IB, IA, TA> self_type;
4276	public:
4277	#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
4278	using matrix_container<self_type>::operator ();
4279	#endif
4280	// ISSUE require type consistency check, is_convertable (IA::size_type, TA::size_type)
4281	typedef typename IA::value_type size_type;
4282	// ISSUE difference_type cannot be deduced for sparse indices, we only know the value_type
4283	typedef std::ptrdiff_t difference_type;
4284	// size_type for the data arrays.
4285	typedef typename IA::size_type array_size_type;
4286	typedef T value_type;
4287	typedef const T &const_reference;
4288	#ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
4289	typedef T &reference;
4290	#else
4291	typedef sparse_matrix_element<self_type> reference;
4292	#endif
4293	typedef IA index_array_type;
4294	typedef TA value_array_type;
4295	typedef const matrix_reference<const self_type> const_closure_type;
4296	typedef matrix_reference<self_type> closure_type;
4297	typedef coordinate_vector<T, IB, IA, TA> vector_temporary_type;
4298	typedef self_type matrix_temporary_type;
4299	typedef sparse_tag storage_category;
4300	typedef typename L::orientation_category orientation_category;
4301
4302	// Construction and destruction
4303	BOOST_UBLAS_INLINE
4304	coordinate_matrix ():
4305	matrix_container<self_type> (),
4306	size1_ (0), size2_ (0), capacity_ (restrict_capacity (non_zeros: 0)),
4307	filled_ (0), sorted_filled_ (filled_), sorted_ (true),
4308	index1_data_ (capacity_), index2_data_ (capacity_), value_data_ (capacity_) {
4309	storage_invariants ();
4310	}
4311	BOOST_UBLAS_INLINE
4312	coordinate_matrix (size_type size1, size_type size2, array_size_type non_zeros = 0):
4313	matrix_container<self_type> (),
4314	size1_ (size1), size2_ (size2), capacity_ (restrict_capacity (non_zeros)),
4315	filled_ (0), sorted_filled_ (filled_), sorted_ (true),
4316	index1_data_ (capacity_), index2_data_ (capacity_), value_data_ (capacity_) {
4317	storage_invariants ();
4318	}
4319	BOOST_UBLAS_INLINE
4320	coordinate_matrix (const coordinate_matrix &m):
4321	matrix_container<self_type> (),
4322	size1_ (m.size1_), size2_ (m.size2_), capacity_ (m.capacity_),
4323	filled_ (m.filled_), sorted_filled_ (m.sorted_filled_), sorted_ (m.sorted_),
4324	index1_data_ (m.index1_data_), index2_data_ (m.index2_data_), value_data_ (m.value_data_) {
4325	storage_invariants ();
4326	}
4327	template<class AE>
4328	BOOST_UBLAS_INLINE
4329	coordinate_matrix (const matrix_expression<AE> &ae, array_size_type non_zeros = 0):
4330	matrix_container<self_type> (),
4331	size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), capacity_ (restrict_capacity (non_zeros)),
4332	filled_ (0), sorted_filled_ (filled_), sorted_ (true),
4333	index1_data_ (capacity_), index2_data_ (capacity_), value_data_ (capacity_) {
4334	storage_invariants ();
4335	matrix_assign<scalar_assign> (*this, ae);
4336	}
4337
4338	// Accessors
4339	BOOST_UBLAS_INLINE
4340	size_type size1 () const {
4341	return size1_;
4342	}
4343	BOOST_UBLAS_INLINE
4344	size_type size2 () const {
4345	return size2_;
4346	}
4347	BOOST_UBLAS_INLINE
4348	size_type nnz_capacity () const {
4349	return capacity_;
4350	}
4351	BOOST_UBLAS_INLINE
4352	size_type nnz () const {
4353	return filled_;
4354	}
4355
4356	// Storage accessors
4357	BOOST_UBLAS_INLINE
4358	static size_type index_base () {
4359	return IB;
4360	}
4361	BOOST_UBLAS_INLINE
4362	array_size_type filled () const {
4363	return filled_;
4364	}
4365	BOOST_UBLAS_INLINE
4366	const index_array_type &index1_data () const {
4367	return index1_data_;
4368	}
4369	BOOST_UBLAS_INLINE
4370	const index_array_type &index2_data () const {
4371	return index2_data_;
4372	}
4373	BOOST_UBLAS_INLINE
4374	const value_array_type &value_data () const {
4375	return value_data_;
4376	}
4377	BOOST_UBLAS_INLINE
4378	void set_filled (const array_size_type &filled) {
4379	// Make sure that storage_invariants() succeeds
4380	if (sorted_ && filled < filled_)
4381	sorted_filled_ = filled;
4382	else
4383	sorted_ = (sorted_filled_ == filled);
4384	filled_ = filled;
4385	storage_invariants ();
4386	}
4387	BOOST_UBLAS_INLINE
4388	index_array_type &index1_data () {
4389	return index1_data_;
4390	}
4391	BOOST_UBLAS_INLINE
4392	index_array_type &index2_data () {
4393	return index2_data_;
4394	}
4395	BOOST_UBLAS_INLINE
4396	value_array_type &value_data () {
4397	return value_data_;
4398	}
4399
4400	// Resizing
4401	private:
4402	BOOST_UBLAS_INLINE
4403	array_size_type restrict_capacity (array_size_type non_zeros) const {
4404	// minimum non_zeros
4405	non_zeros = (std::max) (non_zeros, array_size_type((std::min) (size1_, size2_)));
4406	// ISSUE no maximum as coordinate may contain inserted duplicates
4407	return non_zeros;
4408	}
4409	public:
4410	BOOST_UBLAS_INLINE
4411	void resize (size_type size1, size_type size2, bool preserve = true) {
4412	// FIXME preserve unimplemented
4413	BOOST_UBLAS_CHECK (!preserve, internal_logic ());
4414	size1_ = size1;
4415	size2_ = size2;
4416	capacity_ = restrict_capacity (non_zeros: capacity_);
4417	index1_data_.resize (capacity_);
4418	index2_data_.resize (capacity_);
4419	value_data_.resize (capacity_);
4420	filled_ = 0;
4421	sorted_filled_ = filled_;
4422	sorted_ = true;
4423	storage_invariants ();
4424	}
4425
4426	// Reserving
4427	BOOST_UBLAS_INLINE
4428	void reserve (array_size_type non_zeros, bool preserve = true) {
4429	sort (); // remove duplicate elements
4430	capacity_ = restrict_capacity (non_zeros);
4431	if (preserve) {
4432	index1_data_.resize (capacity_, size_type ());
4433	index2_data_.resize (capacity_, size_type ());
4434	value_data_.resize (capacity_, value_type ());
4435	filled_ = (std::min) (capacity_, filled_);
4436	}
4437	else {
4438	index1_data_.resize (capacity_);
4439	index2_data_.resize (capacity_);
4440	value_data_.resize (capacity_);
4441	filled_ = 0;
4442	}
4443	sorted_filled_ = filled_;
4444	storage_invariants ();
4445	}
4446
4447	// Element support
4448	BOOST_UBLAS_INLINE
4449	pointer find_element (size_type i, size_type j) {
4450	return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i, j));
4451	}
4452	BOOST_UBLAS_INLINE
4453	const_pointer find_element (size_type i, size_type j) const {
4454	sort ();
4455	size_type element1 (layout_type::index_M (i, j));
4456	size_type element2 (layout_type::index_m (i, j));
4457	vector_const_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (zero_based_index: element1), std::less<size_type> ()));
4458	vector_const_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (zero_based_index: element1), std::less<size_type> ()));
4459	if (itv_begin == itv_end)
4460	return 0;
4461	const_subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ()));
4462	const_subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ()));
4463	const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (zero_based_index: element2), std::less<size_type> ()));
4464	if (it == it_end || *it != k_based (zero_based_index: element2))
4465	return 0;
4466	return &value_data_ [it - index2_data_.begin ()];
4467	}
4468
4469	// Element access
4470	BOOST_UBLAS_INLINE
4471	const_reference operator () (size_type i, size_type j) const {
4472	const_pointer p = find_element (i, j);
4473	if (p)
4474	return *p;
4475	else
4476	return zero_;
4477	}
4478	BOOST_UBLAS_INLINE
4479	reference operator () (size_type i, size_type j) {
4480	#ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
4481	pointer p = find_element (i, j);
4482	if (p)
4483	return *p;
4484	else
4485	return insert_element (i, j, value_type/*zero*/());
4486	#else
4487	return reference (*this, i, j);
4488	#endif
4489	}
4490
4491	// Element assignment
4492	BOOST_UBLAS_INLINE
4493	void append_element (size_type i, size_type j, const_reference t) {
4494	if (filled_ >= capacity_)
4495	reserve (non_zeros: 2 * filled_, preserve: true);
4496	BOOST_UBLAS_CHECK (filled_ < capacity_, internal_logic ());
4497	size_type element1 = layout_type::index_M (i, j);
4498	size_type element2 = layout_type::index_m (i, j);
4499	index1_data_ [filled_] = k_based (zero_based_index: element1);
4500	index2_data_ [filled_] = k_based (zero_based_index: element2);
4501	value_data_ [filled_] = t;
4502	++ filled_;
4503	sorted_ = false;
4504	storage_invariants ();
4505	}
4506	BOOST_UBLAS_INLINE
4507	true_reference insert_element (size_type i, size_type j, const_reference t) {
4508	BOOST_UBLAS_CHECK (!find_element (i, j), bad_index ()); // duplicate element
4509	append_element (i, j, t);
4510	return value_data_ [filled_ - 1];
4511	}
4512	BOOST_UBLAS_INLINE
4513	void erase_element (size_type i, size_type j) {
4514	size_type element1 = layout_type::index_M (i, j);
4515	size_type element2 = layout_type::index_m (i, j);
4516	sort ();
4517	vector_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (zero_based_index: element1), std::less<size_type> ()));
4518	vector_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (zero_based_index: element1), std::less<size_type> ()));
4519	subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ()));
4520	subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ()));
4521	subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (zero_based_index: element2), std::less<size_type> ()));
4522	if (it != it_end && *it == k_based (zero_based_index: element2)) {
4523	typename std::iterator_traits<subiterator_type>::difference_type n = it - index2_data_.begin ();
4524	vector_subiterator_type itv (index1_data_.begin () + n);
4525	std::copy (itv + 1, index1_data_.begin () + filled_, itv);
4526	std::copy (it + 1, index2_data_.begin () + filled_, it);
4527	typename value_array_type::iterator itt (value_data_.begin () + n);
4528	std::copy (itt + 1, value_data_.begin () + filled_, itt);
4529	-- filled_;
4530	sorted_filled_ = filled_;
4531	}
4532	storage_invariants ();
4533	}
4534
4535	// Zeroing
4536	BOOST_UBLAS_INLINE
4537	void clear () {
4538	filled_ = 0;
4539	sorted_filled_ = filled_;
4540	sorted_ = true;
4541	storage_invariants ();
4542	}
4543
4544	// Assignment
4545	BOOST_UBLAS_INLINE
4546	coordinate_matrix &operator = (const coordinate_matrix &m) {
4547	if (this != &m) {
4548	size1_ = m.size1_;
4549	size2_ = m.size2_;
4550	capacity_ = m.capacity_;
4551	filled_ = m.filled_;
4552	sorted_filled_ = m.sorted_filled_;
4553	sorted_ = m.sorted_;
4554	index1_data_ = m.index1_data_;
4555	index2_data_ = m.index2_data_;
4556	value_data_ = m.value_data_;
4557	BOOST_UBLAS_CHECK (capacity_ == index1_data_.size (), internal_logic ());
4558	BOOST_UBLAS_CHECK (capacity_ == index2_data_.size (), internal_logic ());
4559	BOOST_UBLAS_CHECK (capacity_ == value_data_.size (), internal_logic ());
4560	}
4561	storage_invariants ();
4562	return *this;
4563	}
4564	template<class C> // Container assignment without temporary
4565	BOOST_UBLAS_INLINE
4566	coordinate_matrix &operator = (const matrix_container<C> &m) {
4567	resize (size1: m ().size1 (), size2: m ().size2 (), preserve: false);
4568	assign (m);
4569	return *this;
4570	}
4571	BOOST_UBLAS_INLINE
4572	coordinate_matrix &assign_temporary (coordinate_matrix &m) {
4573	swap (m);
4574	return *this;
4575	}
4576	template<class AE>
4577	BOOST_UBLAS_INLINE
4578	coordinate_matrix &operator = (const matrix_expression<AE> &ae) {
4579	self_type temporary (ae, capacity_);
4580	return assign_temporary (m&: temporary);
4581	}
4582	template<class AE>
4583	BOOST_UBLAS_INLINE
4584	coordinate_matrix &assign (const matrix_expression<AE> &ae) {
4585	matrix_assign<scalar_assign> (*this, ae);
4586	return *this;
4587	}
4588	template<class AE>
4589	BOOST_UBLAS_INLINE
4590	coordinate_matrix& operator += (const matrix_expression<AE> &ae) {
4591	self_type temporary (*this + ae, capacity_);
4592	return assign_temporary (m&: temporary);
4593	}
4594	template<class C> // Container assignment without temporary
4595	BOOST_UBLAS_INLINE
4596	coordinate_matrix &operator += (const matrix_container<C> &m) {
4597	plus_assign (m);
4598	return *this;
4599	}
4600	template<class AE>
4601	BOOST_UBLAS_INLINE
4602	coordinate_matrix &plus_assign (const matrix_expression<AE> &ae) {
4603	matrix_assign<scalar_plus_assign> (*this, ae);
4604	return *this;
4605	}
4606	template<class AE>
4607	BOOST_UBLAS_INLINE
4608	coordinate_matrix& operator -= (const matrix_expression<AE> &ae) {
4609	self_type temporary (*this - ae, capacity_);
4610	return assign_temporary (m&: temporary);
4611	}
4612	template<class C> // Container assignment without temporary
4613	BOOST_UBLAS_INLINE
4614	coordinate_matrix &operator -= (const matrix_container<C> &m) {
4615	minus_assign (m);
4616	return *this;
4617	}
4618	template<class AE>
4619	BOOST_UBLAS_INLINE
4620	coordinate_matrix &minus_assign (const matrix_expression<AE> &ae) {
4621	matrix_assign<scalar_minus_assign> (*this, ae);
4622	return *this;
4623	}
4624	template<class AT>
4625	BOOST_UBLAS_INLINE
4626	coordinate_matrix& operator *= (const AT &at) {
4627	matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
4628	return *this;
4629	}
4630	template<class AT>
4631	BOOST_UBLAS_INLINE
4632	coordinate_matrix& operator /= (const AT &at) {
4633	matrix_assign_scalar<scalar_divides_assign> (*this, at);
4634	return *this;
4635	}
4636
4637	// Swapping
4638	BOOST_UBLAS_INLINE
4639	void swap (coordinate_matrix &m) {
4640	if (this != &m) {
4641	std::swap (size1_, m.size1_);
4642	std::swap (size2_, m.size2_);
4643	std::swap (capacity_, m.capacity_);
4644	std::swap (filled_, m.filled_);
4645	std::swap (sorted_filled_, m.sorted_filled_);
4646	std::swap (sorted_, m.sorted_);
4647	index1_data_.swap (m.index1_data_);
4648	index2_data_.swap (m.index2_data_);
4649	value_data_.swap (m.value_data_);
4650	}
4651	storage_invariants ();
4652	}
4653	BOOST_UBLAS_INLINE
4654	friend void swap (coordinate_matrix &m1, coordinate_matrix &m2) {
4655	m1.swap (m2);
4656	}
4657
4658	// replacement if STL lower bound algorithm for use of inplace_merge
4659	array_size_type lower_bound (array_size_type beg, array_size_type end, array_size_type target) const {
4660	while (end > beg) {
4661	array_size_type mid = (beg + end) / 2;
4662	if (((index1_data_[mid] < index1_data_[target]) ||
4663	((index1_data_[mid] == index1_data_[target]) &&
4664	(index2_data_[mid] < index2_data_[target])))) {
4665	beg = mid + 1;
4666	} else {
4667	end = mid;
4668	}
4669	}
4670	return beg;
4671	}
4672
4673	// specialized replacement of STL inplace_merge to avoid compilation
4674	// problems with respect to the array_triple iterator
4675	void inplace_merge (array_size_type beg, array_size_type mid, array_size_type end) const {
4676	array_size_type len_lef = mid - beg;
4677	array_size_type len_rig = end - mid;
4678
4679	if (len_lef == 1 && len_rig == 1) {
4680	if ((index1_data_[mid] < index1_data_[beg]) ||
4681	((index1_data_[mid] == index1_data_[beg]) && (index2_data_[mid] < index2_data_[beg])))
4682	{
4683	std::swap(index1_data_[beg], index1_data_[mid]);
4684	std::swap(index2_data_[beg], index2_data_[mid]);
4685	std::swap(value_data_[beg], value_data_[mid]);
4686	}
4687	} else if (len_lef > 0 && len_rig > 0) {
4688	array_size_type lef_mid, rig_mid;
4689	if (len_lef >= len_rig) {
4690	lef_mid = (beg + mid) / 2;
4691	rig_mid = lower_bound(beg: mid, end, target: lef_mid);
4692	} else {
4693	rig_mid = (mid + end) / 2;
4694	lef_mid = lower_bound(beg, end: mid, target: rig_mid);
4695	}
4696	std::rotate(&index1_data_[0] + lef_mid, &index1_data_[0] + mid, &index1_data_[0] + rig_mid);
4697	std::rotate(&index2_data_[0] + lef_mid, &index2_data_[0] + mid, &index2_data_[0] + rig_mid);
4698	std::rotate(&value_data_[0] + lef_mid, &value_data_[0] + mid, &value_data_[0] + rig_mid);
4699
4700	array_size_type new_mid = lef_mid + rig_mid - mid;
4701	inplace_merge(beg, mid: lef_mid, end: new_mid);
4702	inplace_merge(beg: new_mid, mid: rig_mid, end);
4703	}
4704	}
4705
4706	// Sorting and summation of duplicates
4707	BOOST_UBLAS_INLINE
4708	void sort () const {
4709	if (! sorted_ && filled_ > 0) {
4710	typedef index_triple_array<index_array_type, index_array_type, value_array_type> array_triple;
4711	array_triple ita (filled_, index1_data_, index2_data_, value_data_);
4712	#ifndef BOOST_UBLAS_COO_ALWAYS_DO_FULL_SORT
4713	const typename array_triple::iterator iunsorted = ita.begin () + sorted_filled_;
4714	// sort new elements and merge
4715	std::sort (iunsorted, ita.end ());
4716	inplace_merge(beg: 0, mid: sorted_filled_, end: filled_);
4717	#else
4718	const typename array_triple::iterator iunsorted = ita.begin ();
4719	std::sort (iunsorted, ita.end ());
4720	#endif
4721	// sum duplicates with += and remove
4722	array_size_type filled = 0;
4723	for (array_size_type i = 1; i < filled_; ++ i) {
4724	if (index1_data_ [filled] != index1_data_ [i] ||
4725	index2_data_ [filled] != index2_data_ [i]) {
4726	++ filled;
4727	if (filled != i) {
4728	index1_data_ [filled] = index1_data_ [i];
4729	index2_data_ [filled] = index2_data_ [i];
4730	value_data_ [filled] = value_data_ [i];
4731	}
4732	} else {
4733	value_data_ [filled] += value_data_ [i];
4734	}
4735	}
4736	filled_ = filled + 1;
4737	sorted_filled_ = filled_;
4738	sorted_ = true;
4739	storage_invariants ();
4740	}
4741	}
4742
4743	// Back element insertion and erasure
4744	BOOST_UBLAS_INLINE
4745	void push_back (size_type i, size_type j, const_reference t) {
4746	size_type element1 = layout_type::index_M (i, j);
4747	size_type element2 = layout_type::index_m (i, j);
4748	// must maintain sort order
4749	BOOST_UBLAS_CHECK (sorted_ &&
4750	(filled_ == 0 ||
4751	index1_data_ [filled_ - 1] < k_based (element1) ||
4752	(index1_data_ [filled_ - 1] == k_based (element1) && index2_data_ [filled_ - 1] < k_based (element2)))
4753	, external_logic ());
4754	if (filled_ >= capacity_)
4755	reserve (non_zeros: 2 * filled_, preserve: true);
4756	BOOST_UBLAS_CHECK (filled_ < capacity_, internal_logic ());
4757	index1_data_ [filled_] = k_based (zero_based_index: element1);
4758	index2_data_ [filled_] = k_based (zero_based_index: element2);
4759	value_data_ [filled_] = t;
4760	++ filled_;
4761	sorted_filled_ = filled_;
4762	storage_invariants ();
4763	}
4764	BOOST_UBLAS_INLINE
4765	void pop_back () {
4766	// ISSUE invariants could be simpilfied if sorted required as precondition
4767	BOOST_UBLAS_CHECK (filled_ > 0, external_logic ());
4768	-- filled_;
4769	sorted_filled_ = (std::min) (sorted_filled_, filled_);
4770	sorted_ = sorted_filled_ = filled_;
4771	storage_invariants ();
4772	}
4773
4774	// Iterator types
4775	private:
4776	// Use index array iterator
4777	typedef typename IA::const_iterator vector_const_subiterator_type;
4778	typedef typename IA::iterator vector_subiterator_type;
4779	typedef typename IA::const_iterator const_subiterator_type;
4780	typedef typename IA::iterator subiterator_type;
4781
4782	BOOST_UBLAS_INLINE
4783	true_reference at_element (size_type i, size_type j) {
4784	pointer p = find_element (i, j);
4785	BOOST_UBLAS_CHECK (p, bad_index ());
4786	return *p;
4787	}
4788
4789	public:
4790	class const_iterator1;
4791	class iterator1;
4792	class const_iterator2;
4793	class iterator2;
4794	typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
4795	typedef reverse_iterator_base1<iterator1> reverse_iterator1;
4796	typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
4797	typedef reverse_iterator_base2<iterator2> reverse_iterator2;
4798
4799	// Element lookup
4800	// BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
4801	const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const {
4802	sort ();
4803	for (;;) {
4804	size_type address1 (layout_type::index_M (i, j));
4805	size_type address2 (layout_type::index_m (i, j));
4806	vector_const_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (zero_based_index: address1), std::less<size_type> ()));
4807	vector_const_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (zero_based_index: address1), std::less<size_type> ()));
4808
4809	const_subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ()));
4810	const_subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ()));
4811
4812	const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (zero_based_index: address2), std::less<size_type> ()));
4813	vector_const_subiterator_type itv (index1_data_.begin () + (it - index2_data_.begin ()));
4814	if (rank == 0)
4815	return const_iterator1 (*this, rank, i, j, itv, it);
4816	if (it != it_end && zero_based (k_based_index: *it) == address2)
4817	return const_iterator1 (*this, rank, i, j, itv, it);
4818	if (direction > 0) {
4819	if (layout_type::fast_i ()) {
4820	if (it == it_end)
4821	return const_iterator1 (*this, rank, i, j, itv, it);
4822	i = zero_based (k_based_index: *it);
4823	} else {
4824	if (i >= size1_)
4825	return const_iterator1 (*this, rank, i, j, itv, it);
4826	++ i;
4827	}
4828	} else /* if (direction < 0) */ {
4829	if (layout_type::fast_i ()) {
4830	if (it == index2_data_.begin () + array_size_type (zero_based (k_based_index: *itv)))
4831	return const_iterator1 (*this, rank, i, j, itv, it);
4832	i = zero_based (k_based_index: *(it - 1));
4833	} else {
4834	if (i == 0)
4835	return const_iterator1 (*this, rank, i, j, itv, it);
4836	-- i;
4837	}
4838	}
4839	}
4840	}
4841	// BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
4842	iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) {
4843	sort ();
4844	for (;;) {
4845	size_type address1 (layout_type::index_M (i, j));
4846	size_type address2 (layout_type::index_m (i, j));
4847	vector_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (zero_based_index: address1), std::less<size_type> ()));
4848	vector_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (zero_based_index: address1), std::less<size_type> ()));
4849
4850	subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ()));
4851	subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ()));
4852
4853	subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (zero_based_index: address2), std::less<size_type> ()));
4854	vector_subiterator_type itv (index1_data_.begin () + (it - index2_data_.begin ()));
4855	if (rank == 0)
4856	return iterator1 (*this, rank, i, j, itv, it);
4857	if (it != it_end && zero_based (k_based_index: *it) == address2)
4858	return iterator1 (*this, rank, i, j, itv, it);
4859	if (direction > 0) {
4860	if (layout_type::fast_i ()) {
4861	if (it == it_end)
4862	return iterator1 (*this, rank, i, j, itv, it);
4863	i = zero_based (k_based_index: *it);
4864	} else {
4865	if (i >= size1_)
4866	return iterator1 (*this, rank, i, j, itv, it);
4867	++ i;
4868	}
4869	} else /* if (direction < 0) */ {
4870	if (layout_type::fast_i ()) {
4871	if (it == index2_data_.begin () + array_size_type (zero_based (k_based_index: *itv)))
4872	return iterator1 (*this, rank, i, j, itv, it);
4873	i = zero_based (k_based_index: *(it - 1));
4874	} else {
4875	if (i == 0)
4876	return iterator1 (*this, rank, i, j, itv, it);
4877	-- i;
4878	}
4879	}
4880	}
4881	}
4882	// BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
4883	const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const {
4884	sort ();
4885	for (;;) {
4886	size_type address1 (layout_type::index_M (i, j));
4887	size_type address2 (layout_type::index_m (i, j));
4888	vector_const_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (zero_based_index: address1), std::less<size_type> ()));
4889	vector_const_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (zero_based_index: address1), std::less<size_type> ()));
4890
4891	const_subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ()));
4892	const_subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ()));
4893
4894	const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (zero_based_index: address2), std::less<size_type> ()));
4895	vector_const_subiterator_type itv (index1_data_.begin () + (it - index2_data_.begin ()));
4896	if (rank == 0)
4897	return const_iterator2 (*this, rank, i, j, itv, it);
4898	if (it != it_end && zero_based (k_based_index: *it) == address2)
4899	return const_iterator2 (*this, rank, i, j, itv, it);
4900	if (direction > 0) {
4901	if (layout_type::fast_j ()) {
4902	if (it == it_end)
4903	return const_iterator2 (*this, rank, i, j, itv, it);
4904	j = zero_based (k_based_index: *it);
4905	} else {
4906	if (j >= size2_)
4907	return const_iterator2 (*this, rank, i, j, itv, it);
4908	++ j;
4909	}
4910	} else /* if (direction < 0) */ {
4911	if (layout_type::fast_j ()) {
4912	if (it == index2_data_.begin () + array_size_type (zero_based (k_based_index: *itv)))
4913	return const_iterator2 (*this, rank, i, j, itv, it);
4914	j = zero_based (k_based_index: *(it - 1));
4915	} else {
4916	if (j == 0)
4917	return const_iterator2 (*this, rank, i, j, itv, it);
4918	-- j;
4919	}
4920	}
4921	}
4922	}
4923	// BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
4924	iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) {
4925	sort ();
4926	for (;;) {
4927	size_type address1 (layout_type::index_M (i, j));
4928	size_type address2 (layout_type::index_m (i, j));
4929	vector_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (zero_based_index: address1), std::less<size_type> ()));
4930	vector_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (zero_based_index: address1), std::less<size_type> ()));
4931
4932	subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ()));
4933	subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ()));
4934
4935	subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (zero_based_index: address2), std::less<size_type> ()));
4936	vector_subiterator_type itv (index1_data_.begin () + (it - index2_data_.begin ()));
4937	if (rank == 0)
4938	return iterator2 (*this, rank, i, j, itv, it);
4939	if (it != it_end && zero_based (k_based_index: *it) == address2)
4940	return iterator2 (*this, rank, i, j, itv, it);
4941	if (direction > 0) {
4942	if (layout_type::fast_j ()) {
4943	if (it == it_end)
4944	return iterator2 (*this, rank, i, j, itv, it);
4945	j = zero_based (k_based_index: *it);
4946	} else {
4947	if (j >= size2_)
4948	return iterator2 (*this, rank, i, j, itv, it);
4949	++ j;
4950	}
4951	} else /* if (direction < 0) */ {
4952	if (layout_type::fast_j ()) {
4953	if (it == index2_data_.begin () + array_size_type (zero_based (k_based_index: *itv)))
4954	return iterator2 (*this, rank, i, j, itv, it);
4955	j = zero_based (k_based_index: *(it - 1));
4956	} else {
4957	if (j == 0)
4958	return iterator2 (*this, rank, i, j, itv, it);
4959	-- j;
4960	}
4961	}
4962	}
4963	}
4964
4965
4966	class const_iterator1:
4967	public container_const_reference<coordinate_matrix>,
4968	public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
4969	const_iterator1, value_type> {
4970	public:
4971	typedef typename coordinate_matrix::value_type value_type;
4972	typedef typename coordinate_matrix::difference_type difference_type;
4973	typedef typename coordinate_matrix::const_reference reference;
4974	typedef const typename coordinate_matrix::pointer pointer;
4975
4976	typedef const_iterator2 dual_iterator_type;
4977	typedef const_reverse_iterator2 dual_reverse_iterator_type;
4978
4979	// Construction and destruction
4980	BOOST_UBLAS_INLINE
4981	const_iterator1 ():
4982	container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
4983	BOOST_UBLAS_INLINE
4984	const_iterator1 (const self_type &m, int rank, size_type i, size_type j, const vector_const_subiterator_type &itv, const const_subiterator_type &it):
4985	container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
4986	BOOST_UBLAS_INLINE
4987	const_iterator1 (const iterator1 &it):
4988	container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {}
4989
4990	// Arithmetic
4991	BOOST_UBLAS_INLINE
4992	const_iterator1 &operator ++ () {
4993	if (rank_ == 1 && layout_type::fast_i ())
4994	++ it_;
4995	else {
4996	i_ = index1 () + 1;
4997	if (rank_ == 1)
4998	*this = (*this) ().find1 (rank_, i_, j_, 1);
4999	}
5000	return *this;
5001	}
5002	BOOST_UBLAS_INLINE
5003	const_iterator1 &operator -- () {
5004	if (rank_ == 1 && layout_type::fast_i ())
5005	-- it_;
5006	else {
5007	i_ = index1 () - 1;
5008	if (rank_ == 1)
5009	*this = (*this) ().find1 (rank_, i_, j_, -1);
5010	}
5011	return *this;
5012	}
5013
5014	// Dereference
5015	BOOST_UBLAS_INLINE
5016	const_reference operator * () const {
5017	BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
5018	BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
5019	if (rank_ == 1) {
5020	return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
5021	} else {
5022	return (*this) () (i_, j_);
5023	}
5024	}
5025
5026	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
5027	BOOST_UBLAS_INLINE
5028	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5029	typename self_type::
5030	#endif
5031	const_iterator2 begin () const {
5032	const self_type &m = (*this) ();
5033	return m.find2 (1, index1 (), 0);
5034	}
5035	BOOST_UBLAS_INLINE
5036	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5037	typename self_type::
5038	#endif
5039	const_iterator2 cbegin () const {
5040	return begin ();
5041	}
5042	BOOST_UBLAS_INLINE
5043	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5044	typename self_type::
5045	#endif
5046	const_iterator2 end () const {
5047	const self_type &m = (*this) ();
5048	return m.find2 (1, index1 (), m.size2 ());
5049	}
5050	BOOST_UBLAS_INLINE
5051	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5052	typename self_type::
5053	#endif
5054	const_iterator2 cend () const {
5055	return end ();
5056	}
5057	BOOST_UBLAS_INLINE
5058	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5059	typename self_type::
5060	#endif
5061	const_reverse_iterator2 rbegin () const {
5062	return const_reverse_iterator2 (end ());
5063	}
5064	BOOST_UBLAS_INLINE
5065	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5066	typename self_type::
5067	#endif
5068	const_reverse_iterator2 crbegin () const {
5069	return rbegin ();
5070	}
5071	BOOST_UBLAS_INLINE
5072	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5073	typename self_type::
5074	#endif
5075	const_reverse_iterator2 rend () const {
5076	return const_reverse_iterator2 (begin ());
5077	}
5078	BOOST_UBLAS_INLINE
5079	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5080	typename self_type::
5081	#endif
5082	const_reverse_iterator2 crend () const {
5083	return rend ();
5084	}
5085	#endif
5086
5087	// Indices
5088	BOOST_UBLAS_INLINE
5089	size_type index1 () const {
5090	BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
5091	if (rank_ == 1) {
5092	BOOST_UBLAS_CHECK (layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
5093	return layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
5094	} else {
5095	return i_;
5096	}
5097	}
5098	BOOST_UBLAS_INLINE
5099	size_type index2 () const {
5100	if (rank_ == 1) {
5101	BOOST_UBLAS_CHECK (layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
5102	return layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
5103	} else {
5104	return j_;
5105	}
5106	}
5107
5108	// Assignment
5109	BOOST_UBLAS_INLINE
5110	const_iterator1 &operator = (const const_iterator1 &it) {
5111	container_const_reference<self_type>::assign (&it ());
5112	rank_ = it.rank_;
5113	i_ = it.i_;
5114	j_ = it.j_;
5115	itv_ = it.itv_;
5116	it_ = it.it_;
5117	return *this;
5118	}
5119
5120	// Comparison
5121	BOOST_UBLAS_INLINE
5122	bool operator == (const const_iterator1 &it) const {
5123	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
5124	// BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
5125	if (rank_ == 1 || it.rank_ == 1) {
5126	return it_ == it.it_;
5127	} else {
5128	return i_ == it.i_ && j_ == it.j_;
5129	}
5130	}
5131
5132	private:
5133	int rank_;
5134	size_type i_;
5135	size_type j_;
5136	vector_const_subiterator_type itv_;
5137	const_subiterator_type it_;
5138	};
5139
5140	BOOST_UBLAS_INLINE
5141	const_iterator1 begin1 () const {
5142	return find1 (0, 0, 0);
5143	}
5144	BOOST_UBLAS_INLINE
5145	const_iterator1 cbegin1 () const {
5146	return begin1 ();
5147	}
5148	BOOST_UBLAS_INLINE
5149	const_iterator1 end1 () const {
5150	return find1 (0, size1_, 0);
5151	}
5152	BOOST_UBLAS_INLINE
5153	const_iterator1 cend1 () const {
5154	return end1 ();
5155	}
5156
5157	class iterator1:
5158	public container_reference<coordinate_matrix>,
5159	public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
5160	iterator1, value_type> {
5161	public:
5162	typedef typename coordinate_matrix::value_type value_type;
5163	typedef typename coordinate_matrix::difference_type difference_type;
5164	typedef typename coordinate_matrix::true_reference reference;
5165	typedef typename coordinate_matrix::pointer pointer;
5166
5167	typedef iterator2 dual_iterator_type;
5168	typedef reverse_iterator2 dual_reverse_iterator_type;
5169
5170	// Construction and destruction
5171	BOOST_UBLAS_INLINE
5172	iterator1 ():
5173	container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
5174	BOOST_UBLAS_INLINE
5175	iterator1 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it):
5176	container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
5177
5178	// Arithmetic
5179	BOOST_UBLAS_INLINE
5180	iterator1 &operator ++ () {
5181	if (rank_ == 1 && layout_type::fast_i ())
5182	++ it_;
5183	else {
5184	i_ = index1 () + 1;
5185	if (rank_ == 1)
5186	*this = (*this) ().find1 (rank_, i_, j_, 1);
5187	}
5188	return *this;
5189	}
5190	BOOST_UBLAS_INLINE
5191	iterator1 &operator -- () {
5192	if (rank_ == 1 && layout_type::fast_i ())
5193	-- it_;
5194	else {
5195	i_ = index1 () - 1;
5196	if (rank_ == 1)
5197	*this = (*this) ().find1 (rank_, i_, j_, -1);
5198	}
5199	return *this;
5200	}
5201
5202	// Dereference
5203	BOOST_UBLAS_INLINE
5204	reference operator * () const {
5205	BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
5206	BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
5207	if (rank_ == 1) {
5208	return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
5209	} else {
5210	return (*this) ().at_element (i_, j_);
5211	}
5212	}
5213
5214	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
5215	BOOST_UBLAS_INLINE
5216	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5217	typename self_type::
5218	#endif
5219	iterator2 begin () const {
5220	self_type &m = (*this) ();
5221	return m.find2 (1, index1 (), 0);
5222	}
5223	BOOST_UBLAS_INLINE
5224	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5225	typename self_type::
5226	#endif
5227	iterator2 end () const {
5228	self_type &m = (*this) ();
5229	return m.find2 (1, index1 (), m.size2 ());
5230	}
5231	BOOST_UBLAS_INLINE
5232	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5233	typename self_type::
5234	#endif
5235	reverse_iterator2 rbegin () const {
5236	return reverse_iterator2 (end ());
5237	}
5238	BOOST_UBLAS_INLINE
5239	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5240	typename self_type::
5241	#endif
5242	reverse_iterator2 rend () const {
5243	return reverse_iterator2 (begin ());
5244	}
5245	#endif
5246
5247	// Indices
5248	BOOST_UBLAS_INLINE
5249	size_type index1 () const {
5250	BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
5251	if (rank_ == 1) {
5252	BOOST_UBLAS_CHECK (layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
5253	return layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
5254	} else {
5255	return i_;
5256	}
5257	}
5258	BOOST_UBLAS_INLINE
5259	size_type index2 () const {
5260	if (rank_ == 1) {
5261	BOOST_UBLAS_CHECK (layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
5262	return layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
5263	} else {
5264	return j_;
5265	}
5266	}
5267
5268	// Assignment
5269	BOOST_UBLAS_INLINE
5270	iterator1 &operator = (const iterator1 &it) {
5271	container_reference<self_type>::assign (&it ());
5272	rank_ = it.rank_;
5273	i_ = it.i_;
5274	j_ = it.j_;
5275	itv_ = it.itv_;
5276	it_ = it.it_;
5277	return *this;
5278	}
5279
5280	// Comparison
5281	BOOST_UBLAS_INLINE
5282	bool operator == (const iterator1 &it) const {
5283	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
5284	// BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
5285	if (rank_ == 1 || it.rank_ == 1) {
5286	return it_ == it.it_;
5287	} else {
5288	return i_ == it.i_ && j_ == it.j_;
5289	}
5290	}
5291
5292	private:
5293	int rank_;
5294	size_type i_;
5295	size_type j_;
5296	vector_subiterator_type itv_;
5297	subiterator_type it_;
5298
5299	friend class const_iterator1;
5300	};
5301
5302	BOOST_UBLAS_INLINE
5303	iterator1 begin1 () {
5304	return find1 (0, 0, 0);
5305	}
5306	BOOST_UBLAS_INLINE
5307	iterator1 end1 () {
5308	return find1 (0, size1_, 0);
5309	}
5310
5311	class const_iterator2:
5312	public container_const_reference<coordinate_matrix>,
5313	public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
5314	const_iterator2, value_type> {
5315	public:
5316	typedef typename coordinate_matrix::value_type value_type;
5317	typedef typename coordinate_matrix::difference_type difference_type;
5318	typedef typename coordinate_matrix::const_reference reference;
5319	typedef const typename coordinate_matrix::pointer pointer;
5320
5321	typedef const_iterator1 dual_iterator_type;
5322	typedef const_reverse_iterator1 dual_reverse_iterator_type;
5323
5324	// Construction and destruction
5325	BOOST_UBLAS_INLINE
5326	const_iterator2 ():
5327	container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
5328	BOOST_UBLAS_INLINE
5329	const_iterator2 (const self_type &m, int rank, size_type i, size_type j, const vector_const_subiterator_type itv, const const_subiterator_type &it):
5330	container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
5331	BOOST_UBLAS_INLINE
5332	const_iterator2 (const iterator2 &it):
5333	container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {}
5334
5335	// Arithmetic
5336	BOOST_UBLAS_INLINE
5337	const_iterator2 &operator ++ () {
5338	if (rank_ == 1 && layout_type::fast_j ())
5339	++ it_;
5340	else {
5341	j_ = index2 () + 1;
5342	if (rank_ == 1)
5343	*this = (*this) ().find2 (rank_, i_, j_, 1);
5344	}
5345	return *this;
5346	}
5347	BOOST_UBLAS_INLINE
5348	const_iterator2 &operator -- () {
5349	if (rank_ == 1 && layout_type::fast_j ())
5350	-- it_;
5351	else {
5352	j_ = index2 () - 1;
5353	if (rank_ == 1)
5354	*this = (*this) ().find2 (rank_, i_, j_, -1);
5355	}
5356	return *this;
5357	}
5358
5359	// Dereference
5360	BOOST_UBLAS_INLINE
5361	const_reference operator * () const {
5362	BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
5363	BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
5364	if (rank_ == 1) {
5365	return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
5366	} else {
5367	return (*this) () (i_, j_);
5368	}
5369	}
5370
5371	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
5372	BOOST_UBLAS_INLINE
5373	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5374	typename self_type::
5375	#endif
5376	const_iterator1 begin () const {
5377	const self_type &m = (*this) ();
5378	return m.find1 (1, 0, index2 ());
5379	}
5380	BOOST_UBLAS_INLINE
5381	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5382	typename self_type::
5383	#endif
5384	const_iterator1 cbegin () const {
5385	return begin ();
5386	}
5387	BOOST_UBLAS_INLINE
5388	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5389	typename self_type::
5390	#endif
5391	const_iterator1 end () const {
5392	const self_type &m = (*this) ();
5393	return m.find1 (1, m.size1 (), index2 ());
5394	}
5395	BOOST_UBLAS_INLINE
5396	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5397	typename self_type::
5398	#endif
5399	const_iterator1 cend () const {
5400	return end ();
5401	}
5402	BOOST_UBLAS_INLINE
5403	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5404	typename self_type::
5405	#endif
5406	const_reverse_iterator1 rbegin () const {
5407	return const_reverse_iterator1 (end ());
5408	}
5409	BOOST_UBLAS_INLINE
5410	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5411	typename self_type::
5412	#endif
5413	const_reverse_iterator1 crbegin () const {
5414	return rbegin ();
5415	}
5416	BOOST_UBLAS_INLINE
5417	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5418	typename self_type::
5419	#endif
5420	const_reverse_iterator1 rend () const {
5421	return const_reverse_iterator1 (begin ());
5422	}
5423	BOOST_UBLAS_INLINE
5424	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5425	typename self_type::
5426	#endif
5427	const_reverse_iterator1 crend () const {
5428	return rend ();
5429	}
5430	#endif
5431
5432	// Indices
5433	BOOST_UBLAS_INLINE
5434	size_type index1 () const {
5435	if (rank_ == 1) {
5436	BOOST_UBLAS_CHECK (layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
5437	return layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
5438	} else {
5439	return i_;
5440	}
5441	}
5442	BOOST_UBLAS_INLINE
5443	size_type index2 () const {
5444	BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
5445	if (rank_ == 1) {
5446	BOOST_UBLAS_CHECK (layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
5447	return layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
5448	} else {
5449	return j_;
5450	}
5451	}
5452
5453	// Assignment
5454	BOOST_UBLAS_INLINE
5455	const_iterator2 &operator = (const const_iterator2 &it) {
5456	container_const_reference<self_type>::assign (&it ());
5457	rank_ = it.rank_;
5458	i_ = it.i_;
5459	j_ = it.j_;
5460	itv_ = it.itv_;
5461	it_ = it.it_;
5462	return *this;
5463	}
5464
5465	// Comparison
5466	BOOST_UBLAS_INLINE
5467	bool operator == (const const_iterator2 &it) const {
5468	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
5469	// BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
5470	if (rank_ == 1 || it.rank_ == 1) {
5471	return it_ == it.it_;
5472	} else {
5473	return i_ == it.i_ && j_ == it.j_;
5474	}
5475	}
5476
5477	private:
5478	int rank_;
5479	size_type i_;
5480	size_type j_;
5481	vector_const_subiterator_type itv_;
5482	const_subiterator_type it_;
5483	};
5484
5485	BOOST_UBLAS_INLINE
5486	const_iterator2 begin2 () const {
5487	return find2 (0, 0, 0);
5488	}
5489	BOOST_UBLAS_INLINE
5490	const_iterator2 cbegin2 () const {
5491	return begin2 ();
5492	}
5493	BOOST_UBLAS_INLINE
5494	const_iterator2 end2 () const {
5495	return find2 (0, 0, size2_);
5496	}
5497	BOOST_UBLAS_INLINE
5498	const_iterator2 cend2 () const {
5499	return end2 ();
5500	}
5501
5502	class iterator2:
5503	public container_reference<coordinate_matrix>,
5504	public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
5505	iterator2, value_type> {
5506	public:
5507	typedef typename coordinate_matrix::value_type value_type;
5508	typedef typename coordinate_matrix::difference_type difference_type;
5509	typedef typename coordinate_matrix::true_reference reference;
5510	typedef typename coordinate_matrix::pointer pointer;
5511
5512	typedef iterator1 dual_iterator_type;
5513	typedef reverse_iterator1 dual_reverse_iterator_type;
5514
5515	// Construction and destruction
5516	BOOST_UBLAS_INLINE
5517	iterator2 ():
5518	container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {}
5519	BOOST_UBLAS_INLINE
5520	iterator2 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it):
5521	container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {}
5522
5523	// Arithmetic
5524	BOOST_UBLAS_INLINE
5525	iterator2 &operator ++ () {
5526	if (rank_ == 1 && layout_type::fast_j ())
5527	++ it_;
5528	else {
5529	j_ = index2 () + 1;
5530	if (rank_ == 1)
5531	*this = (*this) ().find2 (rank_, i_, j_, 1);
5532	}
5533	return *this;
5534	}
5535	BOOST_UBLAS_INLINE
5536	iterator2 &operator -- () {
5537	if (rank_ == 1 && layout_type::fast_j ())
5538	-- it_;
5539	else {
5540	j_ = index2 ();
5541	if (rank_ == 1)
5542	*this = (*this) ().find2 (rank_, i_, j_, -1);
5543	}
5544	return *this;
5545	}
5546
5547	// Dereference
5548	BOOST_UBLAS_INLINE
5549	reference operator * () const {
5550	BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
5551	BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
5552	if (rank_ == 1) {
5553	return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
5554	} else {
5555	return (*this) ().at_element (i_, j_);
5556	}
5557	}
5558
5559	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
5560	BOOST_UBLAS_INLINE
5561	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5562	typename self_type::
5563	#endif
5564	iterator1 begin () const {
5565	self_type &m = (*this) ();
5566	return m.find1 (1, 0, index2 ());
5567	}
5568	BOOST_UBLAS_INLINE
5569	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5570	typename self_type::
5571	#endif
5572	iterator1 end () const {
5573	self_type &m = (*this) ();
5574	return m.find1 (1, m.size1 (), index2 ());
5575	}
5576	BOOST_UBLAS_INLINE
5577	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5578	typename self_type::
5579	#endif
5580	reverse_iterator1 rbegin () const {
5581	return reverse_iterator1 (end ());
5582	}
5583	BOOST_UBLAS_INLINE
5584	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5585	typename self_type::
5586	#endif
5587	reverse_iterator1 rend () const {
5588	return reverse_iterator1 (begin ());
5589	}
5590	#endif
5591
5592	// Indices
5593	BOOST_UBLAS_INLINE
5594	size_type index1 () const {
5595	if (rank_ == 1) {
5596	BOOST_UBLAS_CHECK (layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ());
5597	return layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
5598	} else {
5599	return i_;
5600	}
5601	}
5602	BOOST_UBLAS_INLINE
5603	size_type index2 () const {
5604	BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
5605	if (rank_ == 1) {
5606	BOOST_UBLAS_CHECK (layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ());
5607	return layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_));
5608	} else {
5609	return j_;
5610	}
5611	}
5612
5613	// Assignment
5614	BOOST_UBLAS_INLINE
5615	iterator2 &operator = (const iterator2 &it) {
5616	container_reference<self_type>::assign (&it ());
5617	rank_ = it.rank_;
5618	i_ = it.i_;
5619	j_ = it.j_;
5620	itv_ = it.itv_;
5621	it_ = it.it_;
5622	return *this;
5623	}
5624
5625	// Comparison
5626	BOOST_UBLAS_INLINE
5627	bool operator == (const iterator2 &it) const {
5628	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
5629	// BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ());
5630	if (rank_ == 1 || it.rank_ == 1) {
5631	return it_ == it.it_;
5632	} else {
5633	return i_ == it.i_ && j_ == it.j_;
5634	}
5635	}
5636
5637	private:
5638	int rank_;
5639	size_type i_;
5640	size_type j_;
5641	vector_subiterator_type itv_;
5642	subiterator_type it_;
5643
5644	friend class const_iterator2;
5645	};
5646
5647	BOOST_UBLAS_INLINE
5648	iterator2 begin2 () {
5649	return find2 (0, 0, 0);
5650	}
5651	BOOST_UBLAS_INLINE
5652	iterator2 end2 () {
5653	return find2 (0, 0, size2_);
5654	}
5655
5656	// Reverse iterators
5657
5658	BOOST_UBLAS_INLINE
5659	const_reverse_iterator1 rbegin1 () const {
5660	return const_reverse_iterator1 (end1 ());
5661	}
5662	BOOST_UBLAS_INLINE
5663	const_reverse_iterator1 crbegin1 () const {
5664	return rbegin1 ();
5665	}
5666	BOOST_UBLAS_INLINE
5667	const_reverse_iterator1 rend1 () const {
5668	return const_reverse_iterator1 (begin1 ());
5669	}
5670	BOOST_UBLAS_INLINE
5671	const_reverse_iterator1 crend1 () const {
5672	return rend1 ();
5673	}
5674
5675	BOOST_UBLAS_INLINE
5676	reverse_iterator1 rbegin1 () {
5677	return reverse_iterator1 (end1 ());
5678	}
5679	BOOST_UBLAS_INLINE
5680	reverse_iterator1 rend1 () {
5681	return reverse_iterator1 (begin1 ());
5682	}
5683
5684	BOOST_UBLAS_INLINE
5685	const_reverse_iterator2 rbegin2 () const {
5686	return const_reverse_iterator2 (end2 ());
5687	}
5688	BOOST_UBLAS_INLINE
5689	const_reverse_iterator2 crbegin2 () const {
5690	return rbegin2 ();
5691	}
5692	BOOST_UBLAS_INLINE
5693	const_reverse_iterator2 rend2 () const {
5694	return const_reverse_iterator2 (begin2 ());
5695	}
5696	BOOST_UBLAS_INLINE
5697	const_reverse_iterator2 crend2 () const {
5698	return rend2 ();
5699	}
5700
5701	BOOST_UBLAS_INLINE
5702	reverse_iterator2 rbegin2 () {
5703	return reverse_iterator2 (end2 ());
5704	}
5705	BOOST_UBLAS_INLINE
5706	reverse_iterator2 rend2 () {
5707	return reverse_iterator2 (begin2 ());
5708	}
5709
5710	// Serialization
5711	template<class Archive>
5712	void serialize(Archive & ar, const unsigned int /* file_version */){
5713	serialization::collection_size_type s1 (size1_);
5714	serialization::collection_size_type s2 (size2_);
5715	ar & serialization::make_nvp(n: "size1",v&: s1);
5716	ar & serialization::make_nvp(n: "size2",v&: s2);
5717	if (Archive::is_loading::value) {
5718	size1_ = s1;
5719	size2_ = s2;
5720	}
5721	ar & serialization::make_nvp("capacity", capacity_);
5722	ar & serialization::make_nvp("filled", filled_);
5723	ar & serialization::make_nvp("sorted_filled", sorted_filled_);
5724	ar & serialization::make_nvp(n: "sorted", v&: sorted_);
5725	ar & serialization::make_nvp("index1_data", index1_data_);
5726	ar & serialization::make_nvp("index2_data", index2_data_);
5727	ar & serialization::make_nvp("value_data", value_data_);
5728	storage_invariants();
5729	}
5730
5731	private:
5732	void storage_invariants () const
5733	{
5734	BOOST_UBLAS_CHECK (capacity_ == index1_data_.size (), internal_logic ());
5735	BOOST_UBLAS_CHECK (capacity_ == index2_data_.size (), internal_logic ());
5736	BOOST_UBLAS_CHECK (capacity_ == value_data_.size (), internal_logic ());
5737	BOOST_UBLAS_CHECK (filled_ <= capacity_, internal_logic ());
5738	BOOST_UBLAS_CHECK (sorted_filled_ <= filled_, internal_logic ());
5739	BOOST_UBLAS_CHECK (sorted_ == (sorted_filled_ == filled_), internal_logic ());
5740	}
5741
5742	size_type size1_;
5743	size_type size2_;
5744	array_size_type capacity_;
5745	mutable array_size_type filled_;
5746	mutable array_size_type sorted_filled_;
5747	mutable bool sorted_;
5748	mutable index_array_type index1_data_;
5749	mutable index_array_type index2_data_;
5750	mutable value_array_type value_data_;
5751	static const value_type zero_;
5752
5753	BOOST_UBLAS_INLINE
5754	static size_type zero_based (size_type k_based_index) {
5755	return k_based_index - IB;
5756	}
5757	BOOST_UBLAS_INLINE
5758	static size_type k_based (size_type zero_based_index) {
5759	return zero_based_index + IB;
5760	}
5761
5762	friend class iterator1;
5763	friend class iterator2;
5764	friend class const_iterator1;
5765	friend class const_iterator2;
5766	};
5767
5768	template<class T, class L, std::size_t IB, class IA, class TA>
5769	const typename coordinate_matrix<T, L, IB, IA, TA>::value_type coordinate_matrix<T, L, IB, IA, TA>::zero_ = value_type/*zero*/();
5770
5771	}}}
5772
5773	#endif
5774