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