Planet

navi

home

PPS

about

screenshots

download

development

forum

Context Navigation

source: downloads/boost_1_34_1/boost/numeric/ublas/matrix_sparse.hpp @ 30

Last change on this file since 30 was 29, checked in by landauf, 17 years ago
updated boost from 1_33_1 to 1_34_1
File size: 206.8 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: