Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

matrix_sparse.hpp @ 857

Revision 857, 211.8 KB checked in by igarcia, 18 years ago (diff)

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

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

Context Navigation

source: NonGTP/Boost/boost/numeric/ublas/matrix_sparse.hpp @ 857

Download in other formats: