1	//
2	// Copyright (c) 2000-2013
3	// Joerg Walter, Mathias Koch, Athanasios Iliopoulos
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_BANDED_
14	#define _BOOST_UBLAS_BANDED_
15
16	#include <boost/numeric/ublas/matrix.hpp>
17	#include <boost/numeric/ublas/detail/temporary.hpp>
18
19	// Iterators based on ideas of Jeremy Siek
20
21	namespace boost { namespace numeric { namespace ublas {
22
23
24	namespace hidden {
25
26
27
28	/** \brief A helper for band_matrix indexing.
29	*
30	* The indexing happens as per the netlib description: http://www.netlib.org/lapack/lug/node124.html.
31	* In the case of a row_major matrix a different approach is followed;
32	*/
33	template <class LayoutType>
34	class banded_indexing { };
35
36	/** \brief A helper for indexing column major banded matrices.
37	*
38	*/
39	template <>
40	class banded_indexing<column_major_tag> {
41	public:
42
43	template <class T>
44	BOOST_UBLAS_INLINE static T size(T /*size1*/, T size2) {
45	return size2;
46	}
47
48	// template <class T>
49	// BOOST_UBLAS_INLINE static bool valid_index(T size1, T /*size2*/, T lower, T upper, T i, T j) {
50	// return (upper+i >= j) && i <= std::min(size1 - 1, j + lower); // upper + i is used by get_index. Maybe find a way to consolidate the operations to increase performance
51	// }
52
53	template <class T>
54	BOOST_UBLAS_INLINE static T get_index(T /*size1*/, T size2, T lower, T upper, T i, T j) {
55	return column_major::element (i: upper + i - j, size_i: lower + 1 + upper, j, size_j: size2);
56	}
57	};
58
59	/** \brief A helper for indexing row major banded matrices.
60	*
61	*/
62	template <>
63	class banded_indexing<row_major_tag> {
64	public:
65
66	template <class T>
67	BOOST_UBLAS_INLINE static T size(T size1, T /*size2*/) {
68	return size1;
69	}
70
71	// template <class T>
72	// BOOST_UBLAS_INLINE static bool valid_index(T /*size1*/, T size2, T lower, T upper, T i, T j) {
73	// return (lower+j >= i) && j <= std::min(size2 - 1, i + upper); // lower + j is used by get_index. Maybe find a way to consolidate the operations to increase performance
74	// }
75
76	template <class T>
77	BOOST_UBLAS_INLINE static T get_index(T size1, T /*size2*/, T lower, T upper, T i, T j) {
78	return row_major::element (i, size_i: size1, j: lower + j - i, size_j: lower + 1 + upper);
79	}
80	};
81
82	}
83
84	/** \brief A banded matrix of values of type \c T.
85	*
86	* For a \f$(mxn)\f$-dimensional banded matrix with \f$l\f$ lower and \f$u\f$ upper diagonals and
87	* \f$0 \leq i < m\f$ and \f$0 \leq j < n\f$, if \f$i>j+l\f$ or \f$i<j-u\f$ then \f$b_{i,j}=0\f$.
88	* The default storage for banded matrices is packed. Orientation and storage can also be specified.
89	* Default is \c row_major and and unbounded_array. It is \b not required by the storage to initialize
90	* elements of the matrix.
91	*
92	* \tparam T the type of object stored in the matrix (like double, float, complex, etc...)
93	* \tparam L the storage organization. It can be either \c row_major or \c column_major. Default is \c row_major
94	* \tparam A the type of Storage array. Default is \c unbounded_array
95	*/
96	template<class T, class L, class A>
97	class banded_matrix:
98	public matrix_container<banded_matrix<T, L, A> > {
99
100	typedef T *pointer;
101	typedef L layout_type;
102	typedef banded_matrix<T, L, A> self_type;
103
104
105
106	public:
107	#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
108	using matrix_container<self_type>::operator ();
109	#endif
110	typedef typename A::size_type size_type;
111	typedef typename A::difference_type difference_type;
112	typedef T value_type;
113	typedef const T &const_reference;
114	typedef T &reference;
115	typedef A array_type;
116	typedef const matrix_reference<const self_type> const_closure_type;
117	typedef matrix_reference<self_type> closure_type;
118	typedef vector<T, A> vector_temporary_type;
119	typedef matrix<T, L, A> matrix_temporary_type; // general sub-matrix
120	typedef packed_tag storage_category;
121	typedef typename L::orientation_category orientation_category;
122
123	private:
124	public:
125
126	// Construction and destruction
127	BOOST_UBLAS_INLINE
128	banded_matrix ():
129	matrix_container<self_type> (),
130	size1_ (0), size2_ (0),
131	lower_ (0), upper_ (0), data_ (0) {}
132	BOOST_UBLAS_INLINE
133	banded_matrix (size_type size1, size_type size2, size_type lower = 0, size_type upper = 0):
134	matrix_container<self_type> (),
135	size1_ (size1), size2_ (size2),
136	lower_ (lower), upper_ (upper),
137	#if defined(BOOST_UBLAS_OWN_BANDED) || (BOOST_UBLAS_LEGACY_BANDED)
138	data_ ((std::max) (size1, size2) * (lower + 1 + upper))
139	#else
140	data_ ( hidden::banded_indexing<orientation_category>::size(size1, size2) * (lower + 1 + upper)) // This is the netlib layout as described here: http://www.netlib.org/lapack/lug/node124.html
141	#endif
142	{
143	}
144	BOOST_UBLAS_INLINE
145	banded_matrix (size_type size1, size_type size2, size_type lower, size_type upper, const array_type &data):
146	matrix_container<self_type> (),
147	size1_ (size1), size2_ (size2),
148	lower_ (lower), upper_ (upper), data_ (data) {}
149	BOOST_UBLAS_INLINE
150	banded_matrix (const banded_matrix &m):
151	matrix_container<self_type> (),
152	size1_ (m.size1_), size2_ (m.size2_),
153	lower_ (m.lower_), upper_ (m.upper_), data_ (m.data_) {}
154	template<class AE>
155	BOOST_UBLAS_INLINE
156	banded_matrix (const matrix_expression<AE> &ae, size_type lower = 0, size_type upper = 0):
157	matrix_container<self_type> (),
158	size1_ (ae ().size1 ()), size2_ (ae ().size2 ()),
159	lower_ (lower), upper_ (upper),
160	#if defined(BOOST_UBLAS_OWN_BANDED) || (BOOST_UBLAS_LEGACY_BANDED)
161	data_ ((std::max) (size1_, size2_) * (lower_ + 1 + upper_))
162	#else
163	data_ ( hidden::banded_indexing<orientation_category>::size(size1_, size2_) * (lower_ + 1 + upper_)) // This is the netlib layout as described here: http://www.netlib.org/lapack/lug/node124.html
164	#endif
165	{
166	matrix_assign<scalar_assign> (*this, ae);
167	}
168
169	// Accessors
170	BOOST_UBLAS_INLINE
171	size_type size1 () const {
172	return size1_;
173	}
174	BOOST_UBLAS_INLINE
175	size_type size2 () const {
176	return size2_;
177	}
178	BOOST_UBLAS_INLINE
179	size_type lower () const {
180	return lower_;
181	}
182	BOOST_UBLAS_INLINE
183	size_type upper () const {
184	return upper_;
185	}
186
187	// Storage accessors
188	BOOST_UBLAS_INLINE
189	const array_type &data () const {
190	return data_;
191	}
192	BOOST_UBLAS_INLINE
193	array_type &data () {
194	return data_;
195	}
196
197	#if !defined (BOOST_UBLAS_OWN_BANDED)||(BOOST_UBLAS_LEGACY_BANDED)
198	BOOST_UBLAS_INLINE
199	bool is_element_in_band(size_type i, size_type j) const{
200	//return (upper_+i >= j) && i <= std::min(size1() - 1, j + lower_); // We don't need to check if i is outside because it is checked anyway in the accessors.
201	return (upper_+i >= j) && i <= ( j + lower_); // Essentially this band has "infinite" positive dimensions
202	}
203	#endif
204	// Resizing
205	BOOST_UBLAS_INLINE
206	void resize (size_type size1, size_type size2, size_type lower = 0, size_type upper = 0, bool preserve = true) {
207	if (preserve) {
208	self_type temporary (size1, size2, lower, upper);
209	detail::matrix_resize_preserve<layout_type> (*this, temporary);
210	}
211	else {
212	data ().resize ((std::max) (size1, size2) * (lower + 1 + upper));
213	size1_ = size1;
214	size2_ = size2;
215	lower_ = lower;
216	upper_ = upper;
217	}
218	}
219
220	BOOST_UBLAS_INLINE
221	void resize_packed_preserve (size_type size1, size_type size2, size_type lower = 0, size_type upper = 0) {
222	size1_ = size1;
223	size2_ = size2;
224	lower_ = lower;
225	upper_ = upper;
226	data ().resize ((std::max) (size1, size2) * (lower + 1 + upper), value_type ());
227	}
228
229	// Element access
230	BOOST_UBLAS_INLINE
231	const_reference operator () (size_type i, size_type j) const {
232	BOOST_UBLAS_CHECK (i < size1_, bad_index ());
233	BOOST_UBLAS_CHECK (j < size2_, bad_index ());
234	#ifdef BOOST_UBLAS_OWN_BANDED
235	const size_type k = (std::max) (i, j);
236	const size_type l = lower_ + j - i;
237	if (k < (std::max) (size1_, size2_) && // TODO: probably use BOOST_UBLAS_CHECK here instead of if
238	l < lower_ + 1 + upper_)
239	return data () [layout_type::element (k, (std::max) (size1_, size2_),
240	l, lower_ + 1 + upper_)];
241	#elif BOOST_UBLAS_LEGACY_BANDED // Prior to version: TODO: add version this is actually incorporated in
242	const size_type k = j;
243	const size_type l = upper_ + i - j;
244	if (k < size2_ &&
245	l < lower_ + 1 + upper_)
246	return data () [layout_type::element (k, size2_,
247	l, lower_ + 1 + upper_)];
248	#else // New default
249	// This is the netlib layout as described here: http://www.netlib.org/lapack/lug/node124.html
250	if ( is_element_in_band( i, j) ) {
251	return data () [hidden::banded_indexing<orientation_category>::get_index(size1_, size2_, lower_, upper_, i, j)];
252	}
253	#endif
254	return zero_;
255	}
256
257	BOOST_UBLAS_INLINE
258	reference at_element (size_type i, size_type j) {
259	BOOST_UBLAS_CHECK (i < size1_, bad_index ());
260	BOOST_UBLAS_CHECK (j < size2_, bad_index ());
261	#ifdef BOOST_UBLAS_OWN_BANDED
262	const size_type k = (std::max) (i, j);
263	const size_type l = lower_ + j - i; // TODO: Don't we need an if or BOOST_UBLAS_CHECK HERE?
264	return data () [layout_type::element (k, (std::max) (size1_, size2_),
265	l, lower_ + 1 + upper_)];
266	#elif BOOST_UBLAS_LEGACY_BANDED // Prior to version: TODO: add version this is actually incorporated in
267	const size_type k = j;
268	const size_type l = upper_ + i - j;
269	if (! (k < size2_ &&
270	l < lower_ + 1 + upper_) ) {
271	bad_index ().raise ();
272	// NEVER reached
273	}
274	return data () [layout_type::element (k, size2_,
275	l, lower_ + 1 + upper_)];
276	#else
277	// This is the netlib layout as described here: http://www.netlib.org/lapack/lug/node124.html
278	BOOST_UBLAS_CHECK(is_element_in_band( i, j) , bad_index());
279	return data () [hidden::banded_indexing<orientation_category>::get_index(size1_, size2_, lower_, upper_, i, j)];
280	#endif
281	}
282	BOOST_UBLAS_INLINE
283	reference operator () (size_type i, size_type j) {
284	BOOST_UBLAS_CHECK (i < size1_, bad_index ());
285	BOOST_UBLAS_CHECK (j < size2_, bad_index ());
286	#ifdef BOOST_UBLAS_OWN_BANDED
287	const size_type k = (std::max) (i, j);
288	const size_type l = lower_ + j - i;
289	if (! (k < (std::max) (size1_, size2_) && // TODO: probably use BOOST_UBLAS_CHECK here instead of if
290	l < lower_ + 1 + upper_) ) {
291	bad_index ().raise ();
292	// NEVER reached
293	}
294	return data () [layout_type::element (k, (std::max) (size1_, size2_),
295	l, lower_ + 1 + upper_)];
296	#elif BOOST_UBLAS_LEGACY_BANDED // Prior to version: TODO: add version this is actually incorporated in
297	const size_type k = j;
298	const size_type l = upper_ + i - j;
299	if (! (k < size2_ &&
300	l < lower_ + 1 + upper_) ) {
301	bad_index ().raise ();
302	// NEVER reached
303	}
304	return data () [layout_type::element (k, size2_,
305	l, lower_ + 1 + upper_)];
306	#else
307	// This is the netlib layout as described here: http://www.netlib.org/lapack/lug/node124.html
308	BOOST_UBLAS_CHECK( is_element_in_band( i, j) , bad_index());
309	return data () [hidden::banded_indexing<orientation_category>::get_index(size1_, size2_, lower_, upper_, i, j)];
310	#endif
311
312	}
313
314	// Element assignment
315	BOOST_UBLAS_INLINE
316	reference insert_element (size_type i, size_type j, const_reference t) {
317	return (operator () (i, j) = t);
318	}
319	BOOST_UBLAS_INLINE
320	void erase_element (size_type i, size_type j) {
321	operator () (i, j) = value_type/*zero*/();
322	}
323
324	// Zeroing
325	BOOST_UBLAS_INLINE
326	void clear () {
327	std::fill (data ().begin (), data ().end (), value_type/*zero*/());
328	}
329
330	// Assignment
331	BOOST_UBLAS_INLINE
332	banded_matrix &operator = (const banded_matrix &m) {
333	size1_ = m.size1_;
334	size2_ = m.size2_;
335	lower_ = m.lower_;
336	upper_ = m.upper_;
337	data () = m.data ();
338	return *this;
339	}
340	BOOST_UBLAS_INLINE
341	banded_matrix &assign_temporary (banded_matrix &m) {
342	swap (m);
343	return *this;
344	}
345	template<class AE>
346	BOOST_UBLAS_INLINE
347	banded_matrix &operator = (const matrix_expression<AE> &ae) {
348	self_type temporary (ae, lower_, upper_);
349	return assign_temporary (m&: temporary);
350	}
351	template<class AE>
352	BOOST_UBLAS_INLINE
353	banded_matrix &assign (const matrix_expression<AE> &ae) {
354	matrix_assign<scalar_assign> (*this, ae);
355	return *this;
356	}
357	template<class AE>
358	BOOST_UBLAS_INLINE
359	banded_matrix& operator += (const matrix_expression<AE> &ae) {
360	self_type temporary (*this + ae, lower_, upper_);
361	return assign_temporary (m&: temporary);
362	}
363	template<class AE>
364	BOOST_UBLAS_INLINE
365	banded_matrix &plus_assign (const matrix_expression<AE> &ae) {
366	matrix_assign<scalar_plus_assign> (*this, ae);
367	return *this;
368	}
369	template<class AE>
370	BOOST_UBLAS_INLINE
371	banded_matrix& operator -= (const matrix_expression<AE> &ae) {
372	self_type temporary (*this - ae, lower_, upper_);
373	return assign_temporary (m&: temporary);
374	}
375	template<class AE>
376	BOOST_UBLAS_INLINE
377	banded_matrix &minus_assign (const matrix_expression<AE> &ae) {
378	matrix_assign<scalar_minus_assign> (*this, ae);
379	return *this;
380	}
381	template<class AT>
382	BOOST_UBLAS_INLINE
383	banded_matrix& operator *= (const AT &at) {
384	matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
385	return *this;
386	}
387	template<class AT>
388	BOOST_UBLAS_INLINE
389	banded_matrix& operator /= (const AT &at) {
390	matrix_assign_scalar<scalar_divides_assign> (*this, at);
391	return *this;
392	}
393
394	// Swapping
395	BOOST_UBLAS_INLINE
396	void swap (banded_matrix &m) {
397	if (this != &m) {
398	std::swap (size1_, m.size1_);
399	std::swap (size2_, m.size2_);
400	std::swap (lower_, m.lower_);
401	std::swap (upper_, m.upper_);
402	data ().swap (m.data ());
403	}
404	}
405	BOOST_UBLAS_INLINE
406	friend void swap (banded_matrix &m1, banded_matrix &m2) {
407	m1.swap (m2);
408	}
409
410	// Iterator types
411	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
412	typedef indexed_iterator1<self_type, packed_random_access_iterator_tag> iterator1;
413	typedef indexed_iterator2<self_type, packed_random_access_iterator_tag> iterator2;
414	typedef indexed_const_iterator1<self_type, packed_random_access_iterator_tag> const_iterator1;
415	typedef indexed_const_iterator2<self_type, packed_random_access_iterator_tag> const_iterator2;
416	#else
417	class const_iterator1;
418	class iterator1;
419	class const_iterator2;
420	class iterator2;
421	#endif
422	typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
423	typedef reverse_iterator_base1<iterator1> reverse_iterator1;
424	typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
425	typedef reverse_iterator_base2<iterator2> reverse_iterator2;
426
427	// Element lookup
428	BOOST_UBLAS_INLINE
429	const_iterator1 find1 (int rank, size_type i, size_type j) const {
430	if (rank == 1) {
431	size_type lower_i = (std::max) (difference_type (j - upper_), difference_type (0));
432	i = (std::max) (i, lower_i);
433	size_type upper_i = (std::min) (j + 1 + lower_, size1_);
434	i = (std::min) (i, upper_i);
435	}
436	return const_iterator1 (*this, i, j);
437	}
438	BOOST_UBLAS_INLINE
439	iterator1 find1 (int rank, size_type i, size_type j) {
440	if (rank == 1) {
441	size_type lower_i = (std::max) (difference_type (j - upper_), difference_type (0));
442	i = (std::max) (i, lower_i);
443	size_type upper_i = (std::min) (j + 1 + lower_, size1_);
444	i = (std::min) (i, upper_i);
445	}
446	return iterator1 (*this, i, j);
447	}
448	BOOST_UBLAS_INLINE
449	const_iterator2 find2 (int rank, size_type i, size_type j) const {
450	if (rank == 1) {
451	size_type lower_j = (std::max) (difference_type (i - lower_), difference_type (0));
452	j = (std::max) (j, lower_j);
453	size_type upper_j = (std::min) (i + 1 + upper_, size2_);
454	j = (std::min) (j, upper_j);
455	}
456	return const_iterator2 (*this, i, j);
457	}
458	BOOST_UBLAS_INLINE
459	iterator2 find2 (int rank, size_type i, size_type j) {
460	if (rank == 1) {
461	size_type lower_j = (std::max) (difference_type (i - lower_), difference_type (0));
462	j = (std::max) (j, lower_j);
463	size_type upper_j = (std::min) (i + 1 + upper_, size2_);
464	j = (std::min) (j, upper_j);
465	}
466	return iterator2 (*this, i, j);
467	}
468
469	// Iterators simply are indices.
470
471	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
472	class const_iterator1:
473	public container_const_reference<banded_matrix>,
474	public random_access_iterator_base<packed_random_access_iterator_tag,
475	const_iterator1, value_type> {
476	public:
477	typedef typename banded_matrix::value_type value_type;
478	typedef typename banded_matrix::difference_type difference_type;
479	typedef typename banded_matrix::const_reference reference;
480	typedef const typename banded_matrix::pointer pointer;
481
482	typedef const_iterator2 dual_iterator_type;
483	typedef const_reverse_iterator2 dual_reverse_iterator_type;
484
485	// Construction and destruction
486	BOOST_UBLAS_INLINE
487	const_iterator1 ():
488	container_const_reference<self_type> (), it1_ (), it2_ () {}
489	BOOST_UBLAS_INLINE
490	const_iterator1 (const self_type &m, size_type it1, size_type it2):
491	container_const_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
492	BOOST_UBLAS_INLINE
493	const_iterator1 (const iterator1 &it):
494	container_const_reference<self_type> (it ()), it1_ (it.it1_), it2_ (it.it2_) {}
495
496	// Arithmetic
497	BOOST_UBLAS_INLINE
498	const_iterator1 &operator ++ () {
499	++ it1_;
500	return *this;
501	}
502	BOOST_UBLAS_INLINE
503	const_iterator1 &operator -- () {
504	-- it1_;
505	return *this;
506	}
507	BOOST_UBLAS_INLINE
508	const_iterator1 &operator += (difference_type n) {
509	it1_ += n;
510	return *this;
511	}
512	BOOST_UBLAS_INLINE
513	const_iterator1 &operator -= (difference_type n) {
514	it1_ -= n;
515	return *this;
516	}
517	BOOST_UBLAS_INLINE
518	difference_type operator - (const const_iterator1 &it) const {
519	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
520	BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
521	return it1_ - it.it1_;
522	}
523
524	// Dereference
525	BOOST_UBLAS_INLINE
526	const_reference operator * () const {
527	return (*this) () (it1_, it2_);
528	}
529	BOOST_UBLAS_INLINE
530	const_reference operator [] (difference_type n) const {
531	return *(*this + n);
532	}
533
534	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
535	BOOST_UBLAS_INLINE
536	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
537	typename self_type::
538	#endif
539	const_iterator2 begin () const {
540	return (*this) ().find2 (1, it1_, 0);
541	}
542	BOOST_UBLAS_INLINE
543	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
544	typename self_type::
545	#endif
546	const_iterator2 cbegin () const {
547	return begin ();
548	}
549	BOOST_UBLAS_INLINE
550	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
551	typename self_type::
552	#endif
553	const_iterator2 end () const {
554	return (*this) ().find2 (1, it1_, (*this) ().size2 ());
555	}
556	BOOST_UBLAS_INLINE
557	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
558	typename self_type::
559	#endif
560	const_iterator2 cend () const {
561	return end ();
562	}
563	BOOST_UBLAS_INLINE
564	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
565	typename self_type::
566	#endif
567	const_reverse_iterator2 rbegin () const {
568	return const_reverse_iterator2 (end ());
569	}
570	BOOST_UBLAS_INLINE
571	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
572	typename self_type::
573	#endif
574	const_reverse_iterator2 crbegin () const {
575	return rbegin ();
576	}
577	BOOST_UBLAS_INLINE
578	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
579	typename self_type::
580	#endif
581	const_reverse_iterator2 rend () const {
582	return const_reverse_iterator2 (begin ());
583	}
584	BOOST_UBLAS_INLINE
585	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
586	typename self_type::
587	#endif
588	const_reverse_iterator2 crend () const {
589	return rend ();
590	}
591	#endif
592
593	// Indices
594	BOOST_UBLAS_INLINE
595	size_type index1 () const {
596	return it1_;
597	}
598	BOOST_UBLAS_INLINE
599	size_type index2 () const {
600	return it2_;
601	}
602
603	// Assignment
604	BOOST_UBLAS_INLINE
605	const_iterator1 &operator = (const const_iterator1 &it) {
606	container_const_reference<self_type>::assign (&it ());
607	it1_ = it.it1_;
608	it2_ = it.it2_;
609	return *this;
610	}
611
612	// Comparison
613	BOOST_UBLAS_INLINE
614	bool operator == (const const_iterator1 &it) const {
615	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
616	BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
617	return it1_ == it.it1_;
618	}
619	BOOST_UBLAS_INLINE
620	bool operator < (const const_iterator1 &it) const {
621	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
622	BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
623	return it1_ < it.it1_;
624	}
625
626	private:
627	size_type it1_;
628	size_type it2_;
629	};
630	#endif
631
632	BOOST_UBLAS_INLINE
633	const_iterator1 begin1 () const {
634	return find1 (0, 0, 0);
635	}
636	BOOST_UBLAS_INLINE
637	const_iterator1 cbegin1 () const {
638	return begin1 ();
639	}
640	BOOST_UBLAS_INLINE
641	const_iterator1 end1 () const {
642	return find1 (0, size1_, 0);
643	}
644	BOOST_UBLAS_INLINE
645	const_iterator1 cend1 () const {
646	return end1 ();
647	}
648
649	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
650	class iterator1:
651	public container_reference<banded_matrix>,
652	public random_access_iterator_base<packed_random_access_iterator_tag,
653	iterator1, value_type> {
654	public:
655	typedef typename banded_matrix::value_type value_type;
656	typedef typename banded_matrix::difference_type difference_type;
657	typedef typename banded_matrix::reference reference;
658	typedef typename banded_matrix::pointer pointer;
659
660	typedef iterator2 dual_iterator_type;
661	typedef reverse_iterator2 dual_reverse_iterator_type;
662
663	// Construction and destruction
664	BOOST_UBLAS_INLINE
665	iterator1 ():
666	container_reference<self_type> (), it1_ (), it2_ () {}
667	BOOST_UBLAS_INLINE
668	iterator1 (self_type &m, size_type it1, size_type it2):
669	container_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
670
671	// Arithmetic
672	BOOST_UBLAS_INLINE
673	iterator1 &operator ++ () {
674	++ it1_;
675	return *this;
676	}
677	BOOST_UBLAS_INLINE
678	iterator1 &operator -- () {
679	-- it1_;
680	return *this;
681	}
682	BOOST_UBLAS_INLINE
683	iterator1 &operator += (difference_type n) {
684	it1_ += n;
685	return *this;
686	}
687	BOOST_UBLAS_INLINE
688	iterator1 &operator -= (difference_type n) {
689	it1_ -= n;
690	return *this;
691	}
692	BOOST_UBLAS_INLINE
693	difference_type operator - (const iterator1 &it) const {
694	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
695	BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
696	return it1_ - it.it1_;
697	}
698
699	// Dereference
700	BOOST_UBLAS_INLINE
701	reference operator * () const {
702	return (*this) ().at_element (it1_, it2_);
703	}
704	BOOST_UBLAS_INLINE
705	reference operator [] (difference_type n) const {
706	return *(*this + n);
707	}
708
709	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
710	BOOST_UBLAS_INLINE
711	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
712	typename self_type::
713	#endif
714	iterator2 begin () const {
715	return (*this) ().find2 (1, it1_, 0);
716	}
717	BOOST_UBLAS_INLINE
718	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
719	typename self_type::
720	#endif
721	iterator2 end () const {
722	return (*this) ().find2 (1, it1_, (*this) ().size2 ());
723	}
724
725	BOOST_UBLAS_INLINE
726	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
727	typename self_type::
728	#endif
729	reverse_iterator2 rbegin () const {
730	return reverse_iterator2 (end ());
731	}
732	BOOST_UBLAS_INLINE
733	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
734	typename self_type::
735	#endif
736	reverse_iterator2 rend () const {
737	return reverse_iterator2 (begin ());
738	}
739	#endif
740
741	// Indices
742	BOOST_UBLAS_INLINE
743	size_type index1 () const {
744	return it1_;
745	}
746	BOOST_UBLAS_INLINE
747	size_type index2 () const {
748	return it2_;
749	}
750
751	// Assignment
752	BOOST_UBLAS_INLINE
753	iterator1 &operator = (const iterator1 &it) {
754	container_reference<self_type>::assign (&it ());
755	it1_ = it.it1_;
756	it2_ = it.it2_;
757	return *this;
758	}
759
760	// Comparison
761	BOOST_UBLAS_INLINE
762	bool operator == (const iterator1 &it) const {
763	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
764	BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
765	return it1_ == it.it1_;
766	}
767	BOOST_UBLAS_INLINE
768	bool operator < (const iterator1 &it) const {
769	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
770	BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
771	return it1_ < it.it1_;
772	}
773
774	private:
775	size_type it1_;
776	size_type it2_;
777
778	friend class const_iterator1;
779	};
780	#endif
781
782	BOOST_UBLAS_INLINE
783	iterator1 begin1 () {
784	return find1 (0, 0, 0);
785	}
786	BOOST_UBLAS_INLINE
787	iterator1 end1 () {
788	return find1 (0, size1_, 0);
789	}
790
791	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
792	class const_iterator2:
793	public container_const_reference<banded_matrix>,
794	public random_access_iterator_base<packed_random_access_iterator_tag,
795	const_iterator2, value_type> {
796	public:
797	typedef typename banded_matrix::value_type value_type;
798	typedef typename banded_matrix::difference_type difference_type;
799	typedef typename banded_matrix::const_reference reference;
800	typedef const typename banded_matrix::pointer pointer;
801
802	typedef const_iterator1 dual_iterator_type;
803	typedef const_reverse_iterator1 dual_reverse_iterator_type;
804
805	// Construction and destruction
806	BOOST_UBLAS_INLINE
807	const_iterator2 ():
808	container_const_reference<self_type> (), it1_ (), it2_ () {}
809	BOOST_UBLAS_INLINE
810	const_iterator2 (const self_type &m, size_type it1, size_type it2):
811	container_const_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
812	BOOST_UBLAS_INLINE
813	const_iterator2 (const iterator2 &it):
814	container_const_reference<self_type> (it ()), it1_ (it.it1_), it2_ (it.it2_) {}
815
816	// Arithmetic
817	BOOST_UBLAS_INLINE
818	const_iterator2 &operator ++ () {
819	++ it2_;
820	return *this;
821	}
822	BOOST_UBLAS_INLINE
823	const_iterator2 &operator -- () {
824	-- it2_;
825	return *this;
826	}
827	BOOST_UBLAS_INLINE
828	const_iterator2 &operator += (difference_type n) {
829	it2_ += n;
830	return *this;
831	}
832	BOOST_UBLAS_INLINE
833	const_iterator2 &operator -= (difference_type n) {
834	it2_ -= n;
835	return *this;
836	}
837	BOOST_UBLAS_INLINE
838	difference_type operator - (const const_iterator2 &it) const {
839	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
840	BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
841	return it2_ - it.it2_;
842	}
843
844	// Dereference
845	BOOST_UBLAS_INLINE
846	const_reference operator * () const {
847	return (*this) () (it1_, it2_);
848	}
849	BOOST_UBLAS_INLINE
850	const_reference operator [] (difference_type n) const {
851	return *(*this + n);
852	}
853
854	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
855	BOOST_UBLAS_INLINE
856	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
857	typename self_type::
858	#endif
859	const_iterator1 begin () const {
860	return (*this) ().find1 (1, 0, it2_);
861	}
862	BOOST_UBLAS_INLINE
863	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
864	typename self_type::
865	#endif
866	const_iterator1 cbegin () const {
867	return begin ();
868	}
869	BOOST_UBLAS_INLINE
870	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
871	typename self_type::
872	#endif
873	const_iterator1 end () const {
874	return (*this) ().find1 (1, (*this) ().size1 (), it2_);
875	}
876	BOOST_UBLAS_INLINE
877	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
878	typename self_type::
879	#endif
880	const_iterator1 cend () const {
881	return end();
882	}
883
884	BOOST_UBLAS_INLINE
885	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
886	typename self_type::
887	#endif
888	const_reverse_iterator1 rbegin () const {
889	return const_reverse_iterator1 (end ());
890	}
891	BOOST_UBLAS_INLINE
892	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
893	typename self_type::
894	#endif
895	const_reverse_iterator1 crbegin () const {
896	return rbegin ();
897	}
898	BOOST_UBLAS_INLINE
899	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
900	typename self_type::
901	#endif
902	const_reverse_iterator1 rend () const {
903	return const_reverse_iterator1 (begin ());
904	}
905	BOOST_UBLAS_INLINE
906	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
907	typename self_type::
908	#endif
909	const_reverse_iterator1 crend () const {
910	return rend ();
911	}
912	#endif
913
914	// Indices
915	BOOST_UBLAS_INLINE
916	size_type index1 () const {
917	return it1_;
918	}
919	BOOST_UBLAS_INLINE
920	size_type index2 () const {
921	return it2_;
922	}
923
924	// Assignment
925	BOOST_UBLAS_INLINE
926	const_iterator2 &operator = (const const_iterator2 &it) {
927	container_const_reference<self_type>::assign (&it ());
928	it1_ = it.it1_;
929	it2_ = it.it2_;
930	return *this;
931	}
932
933	// Comparison
934	BOOST_UBLAS_INLINE
935	bool operator == (const const_iterator2 &it) const {
936	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
937	BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
938	return it2_ == it.it2_;
939	}
940	BOOST_UBLAS_INLINE
941	bool operator < (const const_iterator2 &it) const {
942	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
943	BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
944	return it2_ < it.it2_;
945	}
946
947	private:
948	size_type it1_;
949	size_type it2_;
950	};
951	#endif
952
953	BOOST_UBLAS_INLINE
954	const_iterator2 begin2 () const {
955	return find2 (0, 0, 0);
956	}
957	BOOST_UBLAS_INLINE
958	const_iterator2 cbegin2 () const {
959	return begin2 ();
960	}
961	BOOST_UBLAS_INLINE
962	const_iterator2 end2 () const {
963	return find2 (0, 0, size2_);
964	}
965	BOOST_UBLAS_INLINE
966	const_iterator2 cend2 () const {
967	return end2 ();
968	}
969
970	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
971	class iterator2:
972	public container_reference<banded_matrix>,
973	public random_access_iterator_base<packed_random_access_iterator_tag,
974	iterator2, value_type> {
975	public:
976	typedef typename banded_matrix::value_type value_type;
977	typedef typename banded_matrix::difference_type difference_type;
978	typedef typename banded_matrix::reference reference;
979	typedef typename banded_matrix::pointer pointer;
980
981	typedef iterator1 dual_iterator_type;
982	typedef reverse_iterator1 dual_reverse_iterator_type;
983
984	// Construction and destruction
985	BOOST_UBLAS_INLINE
986	iterator2 ():
987	container_reference<self_type> (), it1_ (), it2_ () {}
988	BOOST_UBLAS_INLINE
989	iterator2 (self_type &m, size_type it1, size_type it2):
990	container_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
991
992	// Arithmetic
993	BOOST_UBLAS_INLINE
994	iterator2 &operator ++ () {
995	++ it2_;
996	return *this;
997	}
998	BOOST_UBLAS_INLINE
999	iterator2 &operator -- () {
1000	-- it2_;
1001	return *this;
1002	}
1003	BOOST_UBLAS_INLINE
1004	iterator2 &operator += (difference_type n) {
1005	it2_ += n;
1006	return *this;
1007	}
1008	BOOST_UBLAS_INLINE
1009	iterator2 &operator -= (difference_type n) {
1010	it2_ -= n;
1011	return *this;
1012	}
1013	BOOST_UBLAS_INLINE
1014	difference_type operator - (const iterator2 &it) const {
1015	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1016	BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
1017	return it2_ - it.it2_;
1018	}
1019
1020	// Dereference
1021	BOOST_UBLAS_INLINE
1022	reference operator * () const {
1023	return (*this) ().at_element (it1_, it2_);
1024	}
1025	BOOST_UBLAS_INLINE
1026	reference operator [] (difference_type n) const {
1027	return *(*this + n);
1028	}
1029
1030	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1031	BOOST_UBLAS_INLINE
1032	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1033	typename self_type::
1034	#endif
1035	iterator1 begin () const {
1036	return (*this) ().find1 (1, 0, it2_);
1037	}
1038	BOOST_UBLAS_INLINE
1039	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1040	typename self_type::
1041	#endif
1042	iterator1 end () const {
1043	return (*this) ().find1 (1, (*this) ().size1 (), it2_);
1044	}
1045	BOOST_UBLAS_INLINE
1046	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1047	typename self_type::
1048	#endif
1049	reverse_iterator1 rbegin () const {
1050	return reverse_iterator1 (end ());
1051	}
1052	BOOST_UBLAS_INLINE
1053	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1054	typename self_type::
1055	#endif
1056	reverse_iterator1 rend () const {
1057	return reverse_iterator1 (begin ());
1058	}
1059	#endif
1060
1061	// Indices
1062	BOOST_UBLAS_INLINE
1063	size_type index1 () const {
1064	return it1_;
1065	}
1066	BOOST_UBLAS_INLINE
1067	size_type index2 () const {
1068	return it2_;
1069	}
1070
1071	// Assignment
1072	BOOST_UBLAS_INLINE
1073	iterator2 &operator = (const iterator2 &it) {
1074	container_reference<self_type>::assign (&it ());
1075	it1_ = it.it1_;
1076	it2_ = it.it2_;
1077	return *this;
1078	}
1079
1080	// Comparison
1081	BOOST_UBLAS_INLINE
1082	bool operator == (const iterator2 &it) const {
1083	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1084	BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
1085	return it2_ == it.it2_;
1086	}
1087	BOOST_UBLAS_INLINE
1088	bool operator < (const iterator2 &it) const {
1089	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1090	BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
1091	return it2_ < it.it2_;
1092	}
1093
1094	private:
1095	size_type it1_;
1096	size_type it2_;
1097
1098	friend class const_iterator2;
1099	};
1100	#endif
1101
1102	BOOST_UBLAS_INLINE
1103	iterator2 begin2 () {
1104	return find2 (0, 0, 0);
1105	}
1106	BOOST_UBLAS_INLINE
1107	iterator2 end2 () {
1108	return find2 (0, 0, size2_);
1109	}
1110
1111	// Reverse iterators
1112
1113	BOOST_UBLAS_INLINE
1114	const_reverse_iterator1 rbegin1 () const {
1115	return const_reverse_iterator1 (end1 ());
1116	}
1117	BOOST_UBLAS_INLINE
1118	const_reverse_iterator1 crbegin1 () const {
1119	return rbegin1 ();
1120	}
1121	BOOST_UBLAS_INLINE
1122	const_reverse_iterator1 rend1 () const {
1123	return const_reverse_iterator1 (begin1 ());
1124	}
1125	BOOST_UBLAS_INLINE
1126	const_reverse_iterator1 crend1 () const {
1127	return rend1 ();
1128	}
1129
1130	BOOST_UBLAS_INLINE
1131	reverse_iterator1 rbegin1 () {
1132	return reverse_iterator1 (end1 ());
1133	}
1134	BOOST_UBLAS_INLINE
1135	reverse_iterator1 rend1 () {
1136	return reverse_iterator1 (begin1 ());
1137	}
1138
1139	BOOST_UBLAS_INLINE
1140	const_reverse_iterator2 rbegin2 () const {
1141	return const_reverse_iterator2 (end2 ());
1142	}
1143	BOOST_UBLAS_INLINE
1144	const_reverse_iterator2 crbegin2 () const {
1145	return rbegin2 ();
1146	}
1147	BOOST_UBLAS_INLINE
1148	const_reverse_iterator2 rend2 () const {
1149	return const_reverse_iterator2 (begin2 ());
1150	}
1151	BOOST_UBLAS_INLINE
1152	const_reverse_iterator2 crend2 () const {
1153	return rend2 ();
1154	}
1155
1156	BOOST_UBLAS_INLINE
1157	reverse_iterator2 rbegin2 () {
1158	return reverse_iterator2 (end2 ());
1159	}
1160	BOOST_UBLAS_INLINE
1161	reverse_iterator2 rend2 () {
1162	return reverse_iterator2 (begin2 ());
1163	}
1164
1165	private:
1166	size_type size1_;
1167	size_type size2_;
1168	size_type lower_;
1169	size_type upper_;
1170	array_type data_;
1171	typedef const value_type const_value_type;
1172	static const_value_type zero_;
1173	};
1174
1175	template<class T, class L, class A>
1176	typename banded_matrix<T, L, A>::const_value_type banded_matrix<T, L, A>::zero_ = value_type/*zero*/();
1177
1178
1179	/** \brief A diagonal matrix of values of type \c T, which is a specialization of a banded matrix
1180	*
1181	* For a \f$(m\times m)\f$-dimensional diagonal matrix, \f$0 \leq i < m\f$ and \f$0 \leq j < m\f$,
1182	* if \f$i\neq j\f$ then \f$b_{i,j}=0\f$. The default storage for diagonal matrices is packed.
1183	* Orientation and storage can also be specified. Default is \c row major \c unbounded_array.
1184	*
1185	* As a specialization of a banded matrix, the constructor of the diagonal matrix creates
1186	* a banded matrix with 0 upper and lower diagonals around the main diagonal and the matrix is
1187	* obviously a square matrix. Operations are optimized based on these 2 assumptions. It is
1188	* \b not required by the storage to initialize elements of the matrix.
1189	*
1190	* \tparam T the type of object stored in the matrix (like double, float, complex, etc...)
1191	* \tparam L the storage organization. It can be either \c row_major or \c column_major. Default is \c row_major
1192	* \tparam A the type of Storage array. Default is \c unbounded_array
1193	*/
1194	template<class T, class L, class A>
1195	class diagonal_matrix:
1196	public banded_matrix<T, L, A> {
1197	public:
1198	typedef typename A::size_type size_type;
1199	typedef banded_matrix<T, L, A> matrix_type;
1200	typedef A array_type;
1201
1202	// Construction and destruction
1203	BOOST_UBLAS_INLINE
1204	diagonal_matrix ():
1205	matrix_type () {}
1206	BOOST_UBLAS_INLINE
1207	diagonal_matrix (size_type size):
1208	matrix_type (size, size) {}
1209	BOOST_UBLAS_INLINE
1210	diagonal_matrix (size_type size, const array_type& data):
1211	matrix_type (size, size, 0, 0, data) {}
1212	BOOST_UBLAS_INLINE
1213	diagonal_matrix (size_type size1, size_type size2):
1214	matrix_type (size1, size2) {}
1215	template<class AE>
1216	BOOST_UBLAS_INLINE
1217	diagonal_matrix (const matrix_expression<AE> &ae):
1218	matrix_type (ae) {}
1219	BOOST_UBLAS_INLINE
1220	~diagonal_matrix () {}
1221
1222	// Assignment
1223	BOOST_UBLAS_INLINE
1224	diagonal_matrix &operator = (const diagonal_matrix &m) {
1225	matrix_type::operator = (m);
1226	return *this;
1227	}
1228	template<class AE>
1229	BOOST_UBLAS_INLINE
1230	diagonal_matrix &operator = (const matrix_expression<AE> &ae) {
1231	matrix_type::operator = (ae);
1232	return *this;
1233	}
1234	};
1235
1236	/** \brief A banded matrix adaptator: convert a any matrix into a banded matrix expression
1237	*
1238	* For a \f$(m\times n)\f$-dimensional matrix, the \c banded_adaptor will provide a banded matrix
1239	* with \f$l\f$ lower and \f$u\f$ upper diagonals and \f$0 \leq i < m\f$ and \f$0 \leq j < n\f$,
1240	* if \f$i>j+l\f$ or \f$i<j-u\f$ then \f$b_{i,j}=0\f$.
1241	*
1242	* Storage and location are based on those of the underlying matrix. This is important because
1243	* a \c banded_adaptor does not copy the matrix data to a new place. Therefore, modifying values
1244	* in a \c banded_adaptor matrix will also modify the underlying matrix too.
1245	*
1246	* \tparam M the type of matrix used to generate a banded matrix
1247	*/
1248	template<class M>
1249	class banded_adaptor:
1250	public matrix_expression<banded_adaptor<M> > {
1251
1252	typedef banded_adaptor<M> self_type;
1253
1254	public:
1255	#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
1256	using matrix_expression<self_type>::operator ();
1257	#endif
1258	typedef const M const_matrix_type;
1259	typedef M matrix_type;
1260	typedef typename M::size_type size_type;
1261	typedef typename M::difference_type difference_type;
1262	typedef typename M::value_type value_type;
1263	typedef typename M::const_reference const_reference;
1264	typedef typename boost::mpl::if_<boost::is_const<M>,
1265	typename M::const_reference,
1266	typename M::reference>::type reference;
1267	typedef typename boost::mpl::if_<boost::is_const<M>,
1268	typename M::const_closure_type,
1269	typename M::closure_type>::type matrix_closure_type;
1270	typedef const self_type const_closure_type;
1271	typedef self_type closure_type;
1272	// Replaced by _temporary_traits to avoid type requirements on M
1273	//typedef typename M::vector_temporary_type vector_temporary_type;
1274	//typedef typename M::matrix_temporary_type matrix_temporary_type;
1275	typedef typename storage_restrict_traits<typename M::storage_category,
1276	packed_proxy_tag>::storage_category storage_category;
1277	typedef typename M::orientation_category orientation_category;
1278
1279	// Construction and destruction
1280	BOOST_UBLAS_INLINE
1281	banded_adaptor (matrix_type &data, size_type lower = 0, size_type upper = 0):
1282	matrix_expression<self_type> (),
1283	data_ (data), lower_ (lower), upper_ (upper) {}
1284	BOOST_UBLAS_INLINE
1285	banded_adaptor (const banded_adaptor &m):
1286	matrix_expression<self_type> (),
1287	data_ (m.data_), lower_ (m.lower_), upper_ (m.upper_) {}
1288
1289	// Accessors
1290	BOOST_UBLAS_INLINE
1291	size_type size1 () const {
1292	return data_.size1 ();
1293	}
1294	BOOST_UBLAS_INLINE
1295	size_type size2 () const {
1296	return data_.size2 ();
1297	}
1298	BOOST_UBLAS_INLINE
1299	size_type lower () const {
1300	return lower_;
1301	}
1302	BOOST_UBLAS_INLINE
1303	size_type upper () const {
1304	return upper_;
1305	}
1306
1307	// Storage accessors
1308	BOOST_UBLAS_INLINE
1309	const matrix_closure_type &data () const {
1310	return data_;
1311	}
1312	BOOST_UBLAS_INLINE
1313	matrix_closure_type &data () {
1314	return data_;
1315	}
1316
1317	#if !defined (BOOST_UBLAS_OWN_BANDED)||(BOOST_UBLAS_LEGACY_BANDED)
1318	BOOST_UBLAS_INLINE
1319	bool is_element_in_band(size_type i, size_type j) const{
1320	//return (upper_+i >= j) && i <= std::min(size1() - 1, j + lower_); // We don't need to check if i is outside because it is checked anyway in the accessors.
1321	return (upper_+i >= j) && i <= ( j + lower_); // Essentially this band has "infinite" positive dimensions
1322	}
1323	#endif
1324
1325	// Element access
1326	#ifndef BOOST_UBLAS_PROXY_CONST_MEMBER
1327	BOOST_UBLAS_INLINE
1328	const_reference operator () (size_type i, size_type j) const {
1329	BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
1330	BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
1331	#ifdef BOOST_UBLAS_OWN_BANDED
1332	size_type k = (std::max) (i, j);
1333	size_type l = lower_ + j - i;
1334	if (k < (std::max) (size1 (), size2 ()) &&
1335	l < lower_ + 1 + upper_)
1336	return data () (i, j);
1337	#elif BOOST_UBLAS_LEGACY_BANDED
1338	size_type k = j;
1339	size_type l = upper_ + i - j;
1340	if (k < size2 () &&
1341	l < lower_ + 1 + upper_)
1342	return data () (i, j);
1343	#else
1344	if (is_element_in_band( i, j))
1345	return data () (i, j);
1346	#endif
1347	return zero_;
1348	}
1349	BOOST_UBLAS_INLINE
1350	reference operator () (size_type i, size_type j) {
1351	BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
1352	BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
1353	#ifdef BOOST_UBLAS_OWN_BANDED
1354	size_type k = (std::max) (i, j);
1355	size_type l = lower_ + j - i;
1356	if (k < (std::max) (size1 (), size2 ()) &&
1357	l < lower_ + 1 + upper_)
1358	return data () (i, j);
1359	#elif BOOST_UBLAS_LEGACY_BANDED
1360	size_type k = j;
1361	size_type l = upper_ + i - j;
1362	if (k < size2 () &&
1363	l < lower_ + 1 + upper_)
1364	return data () (i, j);
1365	#else
1366	if (is_element_in_band( i, j))
1367	return data () (i, j);
1368	#endif
1369	#ifndef BOOST_UBLAS_REFERENCE_CONST_MEMBER
1370	bad_index ().raise ();
1371	#endif
1372	return const_cast<reference>(zero_);
1373	}
1374	#else
1375	BOOST_UBLAS_INLINE
1376	reference operator () (size_type i, size_type j) const {
1377	BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
1378	BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
1379	#ifdef BOOST_UBLAS_OWN_BANDED
1380	size_type k = (std::max) (i, j);
1381	size_type l = lower_ + j - i;
1382	if (k < (std::max) (size1 (), size2 ()) &&
1383	l < lower_ + 1 + upper_)
1384	return data () (i, j);
1385	#elif BOOST_UBLAS_LEGACY_BANDED
1386	size_type k = j;
1387	size_type l = upper_ + i - j;
1388	if (k < size2 () &&
1389	l < lower_ + 1 + upper_)
1390	return data () (i, j);
1391	#else
1392	if (is_element_in_band( i, j))
1393	return data () (i, j);
1394	#endif
1395	#ifndef BOOST_UBLAS_REFERENCE_CONST_MEMBER
1396	bad_index ().raise ();
1397	#endif
1398	return const_cast<reference>(zero_);
1399	}
1400	#endif
1401
1402	// Assignment
1403	BOOST_UBLAS_INLINE
1404	banded_adaptor &operator = (const banded_adaptor &m) {
1405	matrix_assign<scalar_assign> (*this, m);
1406	return *this;
1407	}
1408	BOOST_UBLAS_INLINE
1409	banded_adaptor &assign_temporary (banded_adaptor &m) {
1410	*this = m;
1411	return *this;
1412	}
1413	template<class AE>
1414	BOOST_UBLAS_INLINE
1415	banded_adaptor &operator = (const matrix_expression<AE> &ae) {
1416	matrix_assign<scalar_assign> (*this, matrix<value_type> (ae));
1417	return *this;
1418	}
1419	template<class AE>
1420	BOOST_UBLAS_INLINE
1421	banded_adaptor &assign (const matrix_expression<AE> &ae) {
1422	matrix_assign<scalar_assign> (*this, ae);
1423	return *this;
1424	}
1425	template<class AE>
1426	BOOST_UBLAS_INLINE
1427	banded_adaptor& operator += (const matrix_expression<AE> &ae) {
1428	matrix_assign<scalar_assign> (*this, matrix<value_type> (*this + ae));
1429	return *this;
1430	}
1431	template<class AE>
1432	BOOST_UBLAS_INLINE
1433	banded_adaptor &plus_assign (const matrix_expression<AE> &ae) {
1434	matrix_assign<scalar_plus_assign> (*this, ae);
1435	return *this;
1436	}
1437	template<class AE>
1438	BOOST_UBLAS_INLINE
1439	banded_adaptor& operator -= (const matrix_expression<AE> &ae) {
1440	matrix_assign<scalar_assign> (*this, matrix<value_type> (*this - ae));
1441	return *this;
1442	}
1443	template<class AE>
1444	BOOST_UBLAS_INLINE
1445	banded_adaptor &minus_assign (const matrix_expression<AE> &ae) {
1446	matrix_assign<scalar_minus_assign> (*this, ae);
1447	return *this;
1448	}
1449	template<class AT>
1450	BOOST_UBLAS_INLINE
1451	banded_adaptor& operator *= (const AT &at) {
1452	matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
1453	return *this;
1454	}
1455	template<class AT>
1456	BOOST_UBLAS_INLINE
1457	banded_adaptor& operator /= (const AT &at) {
1458	matrix_assign_scalar<scalar_divides_assign> (*this, at);
1459	return *this;
1460	}
1461
1462	// Closure comparison
1463	BOOST_UBLAS_INLINE
1464	bool same_closure (const banded_adaptor &ba) const {
1465	return (*this).data ().same_closure (ba.data ());
1466	}
1467
1468	// Swapping
1469	BOOST_UBLAS_INLINE
1470	void swap (banded_adaptor &m) {
1471	if (this != &m) {
1472	BOOST_UBLAS_CHECK (lower_ == m.lower_, bad_size ());
1473	BOOST_UBLAS_CHECK (upper_ == m.upper_, bad_size ());
1474	matrix_swap<scalar_swap> (*this, m);
1475	}
1476	}
1477	BOOST_UBLAS_INLINE
1478	friend void swap (banded_adaptor &m1, banded_adaptor &m2) {
1479	m1.swap (m2);
1480	}
1481
1482	// Iterator types
1483	private:
1484	// Use the matrix iterator
1485	typedef typename M::const_iterator1 const_subiterator1_type;
1486	typedef typename boost::mpl::if_<boost::is_const<M>,
1487	typename M::const_iterator1,
1488	typename M::iterator1>::type subiterator1_type;
1489	typedef typename M::const_iterator2 const_subiterator2_type;
1490	typedef typename boost::mpl::if_<boost::is_const<M>,
1491	typename M::const_iterator2,
1492	typename M::iterator2>::type subiterator2_type;
1493
1494	public:
1495	#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1496	typedef indexed_iterator1<self_type, packed_random_access_iterator_tag> iterator1;
1497	typedef indexed_iterator2<self_type, packed_random_access_iterator_tag> iterator2;
1498	typedef indexed_const_iterator1<self_type, packed_random_access_iterator_tag> const_iterator1;
1499	typedef indexed_const_iterator2<self_type, packed_random_access_iterator_tag> const_iterator2;
1500	#else
1501	class const_iterator1;
1502	class iterator1;
1503	class const_iterator2;
1504	class iterator2;
1505	#endif
1506	typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
1507	typedef reverse_iterator_base1<iterator1> reverse_iterator1;
1508	typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
1509	typedef reverse_iterator_base2<iterator2> reverse_iterator2;
1510
1511	// Element lookup
1512	BOOST_UBLAS_INLINE
1513	const_iterator1 find1 (int rank, size_type i, size_type j) const {
1514	if (rank == 1) {
1515	size_type lower_i = (std::max) (difference_type (j - upper_), difference_type (0));
1516	i = (std::max) (i, lower_i);
1517	size_type upper_i = (std::min) (j + 1 + lower_, size1 ());
1518	i = (std::min) (i, upper_i);
1519	}
1520	return const_iterator1 (*this, data ().find1 (rank, i, j));
1521	}
1522	BOOST_UBLAS_INLINE
1523	iterator1 find1 (int rank, size_type i, size_type j) {
1524	if (rank == 1) {
1525	size_type lower_i = (std::max) (difference_type (j - upper_), difference_type (0));
1526	i = (std::max) (i, lower_i);
1527	size_type upper_i = (std::min) (j + 1 + lower_, size1 ());
1528	i = (std::min) (i, upper_i);
1529	}
1530	return iterator1 (*this, data ().find1 (rank, i, j));
1531	}
1532	BOOST_UBLAS_INLINE
1533	const_iterator2 find2 (int rank, size_type i, size_type j) const {
1534	if (rank == 1) {
1535	size_type lower_j = (std::max) (difference_type (i - lower_), difference_type (0));
1536	j = (std::max) (j, lower_j);
1537	size_type upper_j = (std::min) (i + 1 + upper_, size2 ());
1538	j = (std::min) (j, upper_j);
1539	}
1540	return const_iterator2 (*this, data ().find2 (rank, i, j));
1541	}
1542	BOOST_UBLAS_INLINE
1543	iterator2 find2 (int rank, size_type i, size_type j) {
1544	if (rank == 1) {
1545	size_type lower_j = (std::max) (difference_type (i - lower_), difference_type (0));
1546	j = (std::max) (j, lower_j);
1547	size_type upper_j = (std::min) (i + 1 + upper_, size2 ());
1548	j = (std::min) (j, upper_j);
1549	}
1550	return iterator2 (*this, data ().find2 (rank, i, j));
1551	}
1552
1553	// Iterators simply are indices.
1554
1555	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1556	class const_iterator1:
1557	public container_const_reference<banded_adaptor>,
1558	public random_access_iterator_base<typename iterator_restrict_traits<
1559	typename const_subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,
1560	const_iterator1, value_type> {
1561	public:
1562	typedef typename const_subiterator1_type::value_type value_type;
1563	typedef typename const_subiterator1_type::difference_type difference_type;
1564	typedef typename const_subiterator1_type::reference reference;
1565	typedef typename const_subiterator1_type::pointer pointer;
1566
1567	typedef const_iterator2 dual_iterator_type;
1568	typedef const_reverse_iterator2 dual_reverse_iterator_type;
1569
1570	// Construction and destruction
1571	BOOST_UBLAS_INLINE
1572	const_iterator1 ():
1573	container_const_reference<self_type> (), it1_ () {}
1574	BOOST_UBLAS_INLINE
1575	const_iterator1 (const self_type &m, const const_subiterator1_type &it1):
1576	container_const_reference<self_type> (m), it1_ (it1) {}
1577	BOOST_UBLAS_INLINE
1578	const_iterator1 (const iterator1 &it):
1579	container_const_reference<self_type> (it ()), it1_ (it.it1_) {}
1580
1581	// Arithmetic
1582	BOOST_UBLAS_INLINE
1583	const_iterator1 &operator ++ () {
1584	++ it1_;
1585	return *this;
1586	}
1587	BOOST_UBLAS_INLINE
1588	const_iterator1 &operator -- () {
1589	-- it1_;
1590	return *this;
1591	}
1592	BOOST_UBLAS_INLINE
1593	const_iterator1 &operator += (difference_type n) {
1594	it1_ += n;
1595	return *this;
1596	}
1597	BOOST_UBLAS_INLINE
1598	const_iterator1 &operator -= (difference_type n) {
1599	it1_ -= n;
1600	return *this;
1601	}
1602	BOOST_UBLAS_INLINE
1603	difference_type operator - (const const_iterator1 &it) const {
1604	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1605	return it1_ - it.it1_;
1606	}
1607
1608	// Dereference
1609	BOOST_UBLAS_INLINE
1610	const_reference operator * () const {
1611	size_type i = index1 ();
1612	size_type j = index2 ();
1613	BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());
1614	BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());
1615	#ifdef BOOST_UBLAS_OWN_BANDED
1616	size_type k = (std::max) (i, j);
1617	size_type l = (*this) ().lower () + j - i;
1618	if (k < (std::max) ((*this) ().size1 (), (*this) ().size2 ()) &&
1619	l < (*this) ().lower () + 1 + (*this) ().upper ())
1620	return *it1_;
1621	#else
1622	size_type k = j;
1623	size_type l = (*this) ().upper () + i - j;
1624	if (k < (*this) ().size2 () &&
1625	l < (*this) ().lower () + 1 + (*this) ().upper ())
1626	return *it1_;
1627	#endif
1628	return (*this) () (i, j);
1629	}
1630	BOOST_UBLAS_INLINE
1631	const_reference operator [] (difference_type n) const {
1632	return *(*this + n);
1633	}
1634
1635	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1636	BOOST_UBLAS_INLINE
1637	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1638	typename self_type::
1639	#endif
1640	const_iterator2 begin () const {
1641	return (*this) ().find2 (1, index1 (), 0);
1642	}
1643	BOOST_UBLAS_INLINE
1644	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1645	typename self_type::
1646	#endif
1647	const_iterator2 cbegin () const {
1648	return begin ();
1649	}
1650	BOOST_UBLAS_INLINE
1651	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1652	typename self_type::
1653	#endif
1654	const_iterator2 end () const {
1655	return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
1656	}
1657	BOOST_UBLAS_INLINE
1658	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1659	typename self_type::
1660	#endif
1661	const_iterator2 cend () const {
1662	return end ();
1663	}
1664	BOOST_UBLAS_INLINE
1665	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1666	typename self_type::
1667	#endif
1668	const_reverse_iterator2 rbegin () const {
1669	return const_reverse_iterator2 (end ());
1670	}
1671	BOOST_UBLAS_INLINE
1672	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1673	typename self_type::
1674	#endif
1675	const_reverse_iterator2 crbegin () const {
1676	return rbegin ();
1677	}
1678	BOOST_UBLAS_INLINE
1679	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1680	typename self_type::
1681	#endif
1682	const_reverse_iterator2 rend () const {
1683	return const_reverse_iterator2 (begin ());
1684	}
1685	BOOST_UBLAS_INLINE
1686	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1687	typename self_type::
1688	#endif
1689	const_reverse_iterator2 crend () const {
1690	return rend ();
1691	}
1692	#endif
1693
1694	// Indices
1695	BOOST_UBLAS_INLINE
1696	size_type index1 () const {
1697	return it1_.index1 ();
1698	}
1699	BOOST_UBLAS_INLINE
1700	size_type index2 () const {
1701	return it1_.index2 ();
1702	}
1703
1704	// Assignment
1705	BOOST_UBLAS_INLINE
1706	const_iterator1 &operator = (const const_iterator1 &it) {
1707	container_const_reference<self_type>::assign (&it ());
1708	it1_ = it.it1_;
1709	return *this;
1710	}
1711
1712	// Comparison
1713	BOOST_UBLAS_INLINE
1714	bool operator == (const const_iterator1 &it) const {
1715	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1716	return it1_ == it.it1_;
1717	}
1718	BOOST_UBLAS_INLINE
1719	bool operator < (const const_iterator1 &it) const {
1720	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1721	return it1_ < it.it1_;
1722	}
1723
1724	private:
1725	const_subiterator1_type it1_;
1726	};
1727	#endif
1728
1729	BOOST_UBLAS_INLINE
1730	const_iterator1 begin1 () const {
1731	return find1 (0, 0, 0);
1732	}
1733	BOOST_UBLAS_INLINE
1734	const_iterator1 cbegin1 () const {
1735	return begin1 ();
1736	}
1737	BOOST_UBLAS_INLINE
1738	const_iterator1 end1 () const {
1739	return find1 (0, size1 (), 0);
1740	}
1741	BOOST_UBLAS_INLINE
1742	const_iterator1 cend1 () const {
1743	return end1 ();
1744	}
1745
1746	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1747	class iterator1:
1748	public container_reference<banded_adaptor>,
1749	public random_access_iterator_base<typename iterator_restrict_traits<
1750	typename subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,
1751	iterator1, value_type> {
1752	public:
1753	typedef typename subiterator1_type::value_type value_type;
1754	typedef typename subiterator1_type::difference_type difference_type;
1755	typedef typename subiterator1_type::reference reference;
1756	typedef typename subiterator1_type::pointer pointer;
1757
1758	typedef iterator2 dual_iterator_type;
1759	typedef reverse_iterator2 dual_reverse_iterator_type;
1760
1761	// Construction and destruction
1762	BOOST_UBLAS_INLINE
1763	iterator1 ():
1764	container_reference<self_type> (), it1_ () {}
1765	BOOST_UBLAS_INLINE
1766	iterator1 (self_type &m, const subiterator1_type &it1):
1767	container_reference<self_type> (m), it1_ (it1) {}
1768
1769	// Arithmetic
1770	BOOST_UBLAS_INLINE
1771	iterator1 &operator ++ () {
1772	++ it1_;
1773	return *this;
1774	}
1775	BOOST_UBLAS_INLINE
1776	iterator1 &operator -- () {
1777	-- it1_;
1778	return *this;
1779	}
1780	BOOST_UBLAS_INLINE
1781	iterator1 &operator += (difference_type n) {
1782	it1_ += n;
1783	return *this;
1784	}
1785	BOOST_UBLAS_INLINE
1786	iterator1 &operator -= (difference_type n) {
1787	it1_ -= n;
1788	return *this;
1789	}
1790	BOOST_UBLAS_INLINE
1791	difference_type operator - (const iterator1 &it) const {
1792	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1793	return it1_ - it.it1_;
1794	}
1795
1796	// Dereference
1797	BOOST_UBLAS_INLINE
1798	reference operator * () const {
1799	size_type i = index1 ();
1800	size_type j = index2 ();
1801	BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());
1802	BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());
1803	#ifdef BOOST_UBLAS_OWN_BANDED
1804	size_type k = (std::max) (i, j);
1805	size_type l = (*this) ().lower () + j - i;
1806	if (k < (std::max) ((*this) ().size1 (), (*this) ().size2 ()) &&
1807	l < (*this) ().lower () + 1 + (*this) ().upper ())
1808	return *it1_;
1809	#else
1810	size_type k = j;
1811	size_type l = (*this) ().upper () + i - j;
1812	if (k < (*this) ().size2 () &&
1813	l < (*this) ().lower () + 1 + (*this) ().upper ())
1814	return *it1_;
1815	#endif
1816	return (*this) () (i, j);
1817	}
1818	BOOST_UBLAS_INLINE
1819	reference operator [] (difference_type n) const {
1820	return *(*this + n);
1821	}
1822
1823	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1824	BOOST_UBLAS_INLINE
1825	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1826	typename self_type::
1827	#endif
1828	iterator2 begin () const {
1829	return (*this) ().find2 (1, index1 (), 0);
1830	}
1831	BOOST_UBLAS_INLINE
1832	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1833	typename self_type::
1834	#endif
1835	iterator2 end () const {
1836	return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
1837	}
1838	BOOST_UBLAS_INLINE
1839	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1840	typename self_type::
1841	#endif
1842	reverse_iterator2 rbegin () const {
1843	return reverse_iterator2 (end ());
1844	}
1845	BOOST_UBLAS_INLINE
1846	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1847	typename self_type::
1848	#endif
1849	reverse_iterator2 rend () const {
1850	return reverse_iterator2 (begin ());
1851	}
1852	#endif
1853
1854	// Indices
1855	BOOST_UBLAS_INLINE
1856	size_type index1 () const {
1857	return it1_.index1 ();
1858	}
1859	BOOST_UBLAS_INLINE
1860	size_type index2 () const {
1861	return it1_.index2 ();
1862	}
1863
1864	// Assignment
1865	BOOST_UBLAS_INLINE
1866	iterator1 &operator = (const iterator1 &it) {
1867	container_reference<self_type>::assign (&it ());
1868	it1_ = it.it1_;
1869	return *this;
1870	}
1871
1872	// Comparison
1873	BOOST_UBLAS_INLINE
1874	bool operator == (const iterator1 &it) const {
1875	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1876	return it1_ == it.it1_;
1877	}
1878	BOOST_UBLAS_INLINE
1879	bool operator < (const iterator1 &it) const {
1880	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1881	return it1_ < it.it1_;
1882	}
1883
1884	private:
1885	subiterator1_type it1_;
1886
1887	friend class const_iterator1;
1888	};
1889	#endif
1890
1891	BOOST_UBLAS_INLINE
1892	iterator1 begin1 () {
1893	return find1 (0, 0, 0);
1894	}
1895	BOOST_UBLAS_INLINE
1896	iterator1 end1 () {
1897	return find1 (0, size1 (), 0);
1898	}
1899
1900	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1901	class const_iterator2:
1902	public container_const_reference<banded_adaptor>,
1903	public random_access_iterator_base<packed_random_access_iterator_tag,
1904	const_iterator2, value_type> {
1905	public:
1906	typedef typename iterator_restrict_traits<typename const_subiterator2_type::iterator_category,
1907	packed_random_access_iterator_tag>::iterator_category iterator_category;
1908	typedef typename const_subiterator2_type::value_type value_type;
1909	typedef typename const_subiterator2_type::difference_type difference_type;
1910	typedef typename const_subiterator2_type::reference reference;
1911	typedef typename const_subiterator2_type::pointer pointer;
1912
1913	typedef const_iterator1 dual_iterator_type;
1914	typedef const_reverse_iterator1 dual_reverse_iterator_type;
1915
1916	// Construction and destruction
1917	BOOST_UBLAS_INLINE
1918	const_iterator2 ():
1919	container_const_reference<self_type> (), it2_ () {}
1920	BOOST_UBLAS_INLINE
1921	const_iterator2 (const self_type &m, const const_subiterator2_type &it2):
1922	container_const_reference<self_type> (m), it2_ (it2) {}
1923	BOOST_UBLAS_INLINE
1924	const_iterator2 (const iterator2 &it):
1925	container_const_reference<self_type> (it ()), it2_ (it.it2_) {}
1926
1927	// Arithmetic
1928	BOOST_UBLAS_INLINE
1929	const_iterator2 &operator ++ () {
1930	++ it2_;
1931	return *this;
1932	}
1933	BOOST_UBLAS_INLINE
1934	const_iterator2 &operator -- () {
1935	-- it2_;
1936	return *this;
1937	}
1938	BOOST_UBLAS_INLINE
1939	const_iterator2 &operator += (difference_type n) {
1940	it2_ += n;
1941	return *this;
1942	}
1943	BOOST_UBLAS_INLINE
1944	const_iterator2 &operator -= (difference_type n) {
1945	it2_ -= n;
1946	return *this;
1947	}
1948	BOOST_UBLAS_INLINE
1949	difference_type operator - (const const_iterator2 &it) const {
1950	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1951	return it2_ - it.it2_;
1952	}
1953
1954	// Dereference
1955	BOOST_UBLAS_INLINE
1956	const_reference operator * () const {
1957	size_type i = index1 ();
1958	size_type j = index2 ();
1959	BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());
1960	BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());
1961	#ifdef BOOST_UBLAS_OWN_BANDED
1962	size_type k = (std::max) (i, j);
1963	size_type l = (*this) ().lower () + j - i;
1964	if (k < (std::max) ((*this) ().size1 (), (*this) ().size2 ()) &&
1965	l < (*this) ().lower () + 1 + (*this) ().upper ())
1966	return *it2_;
1967	#else
1968	size_type k = j;
1969	size_type l = (*this) ().upper () + i - j;
1970	if (k < (*this) ().size2 () &&
1971	l < (*this) ().lower () + 1 + (*this) ().upper ())
1972	return *it2_;
1973	#endif
1974	return (*this) () (i, j);
1975	}
1976	BOOST_UBLAS_INLINE
1977	const_reference operator [] (difference_type n) const {
1978	return *(*this + n);
1979	}
1980
1981	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1982	BOOST_UBLAS_INLINE
1983	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1984	typename self_type::
1985	#endif
1986	const_iterator1 begin () const {
1987	return (*this) ().find1 (1, 0, index2 ());
1988	}
1989	BOOST_UBLAS_INLINE
1990	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1991	typename self_type::
1992	#endif
1993	const_iterator1 cbegin () const {
1994	return begin ();
1995	}
1996	BOOST_UBLAS_INLINE
1997	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1998	typename self_type::
1999	#endif
2000	const_iterator1 end () const {
2001	return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
2002	}
2003	BOOST_UBLAS_INLINE
2004	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2005	typename self_type::
2006	#endif
2007	const_iterator1 cend () const {
2008	return end ();
2009	}
2010	BOOST_UBLAS_INLINE
2011	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2012	typename self_type::
2013	#endif
2014	const_reverse_iterator1 rbegin () const {
2015	return const_reverse_iterator1 (end ());
2016	}
2017	BOOST_UBLAS_INLINE
2018	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2019	typename self_type::
2020	#endif
2021	const_reverse_iterator1 crbegin () const {
2022	return rbegin ();
2023	}
2024	BOOST_UBLAS_INLINE
2025	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2026	typename self_type::
2027	#endif
2028	const_reverse_iterator1 rend () const {
2029	return const_reverse_iterator1 (begin ());
2030	}
2031	BOOST_UBLAS_INLINE
2032	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2033	typename self_type::
2034	#endif
2035	const_reverse_iterator1 crend () const {
2036	return rend ();
2037	}
2038	#endif
2039
2040	// Indices
2041	BOOST_UBLAS_INLINE
2042	size_type index1 () const {
2043	return it2_.index1 ();
2044	}
2045	BOOST_UBLAS_INLINE
2046	size_type index2 () const {
2047	return it2_.index2 ();
2048	}
2049
2050	// Assignment
2051	BOOST_UBLAS_INLINE
2052	const_iterator2 &operator = (const const_iterator2 &it) {
2053	container_const_reference<self_type>::assign (&it ());
2054	it2_ = it.it2_;
2055	return *this;
2056	}
2057
2058	// Comparison
2059	BOOST_UBLAS_INLINE
2060	bool operator == (const const_iterator2 &it) const {
2061	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2062	return it2_ == it.it2_;
2063	}
2064	BOOST_UBLAS_INLINE
2065	bool operator < (const const_iterator2 &it) const {
2066	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2067	return it2_ < it.it2_;
2068	}
2069
2070	private:
2071	const_subiterator2_type it2_;
2072	};
2073	#endif
2074
2075	BOOST_UBLAS_INLINE
2076	const_iterator2 begin2 () const {
2077	return find2 (0, 0, 0);
2078	}
2079	BOOST_UBLAS_INLINE
2080	const_iterator2 cbegin2 () const {
2081	return begin2 ();
2082	}
2083	BOOST_UBLAS_INLINE
2084	const_iterator2 end2 () const {
2085	return find2 (0, 0, size2 ());
2086	}
2087	BOOST_UBLAS_INLINE
2088	const_iterator2 cend2 () const {
2089	return end2 ();
2090	}
2091
2092	#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
2093	class iterator2:
2094	public container_reference<banded_adaptor>,
2095	public random_access_iterator_base<typename iterator_restrict_traits<
2096	typename subiterator2_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,
2097	iterator2, value_type> {
2098	public:
2099	typedef typename subiterator2_type::value_type value_type;
2100	typedef typename subiterator2_type::difference_type difference_type;
2101	typedef typename subiterator2_type::reference reference;
2102	typedef typename subiterator2_type::pointer pointer;
2103
2104	typedef iterator1 dual_iterator_type;
2105	typedef reverse_iterator1 dual_reverse_iterator_type;
2106
2107	// Construction and destruction
2108	BOOST_UBLAS_INLINE
2109	iterator2 ():
2110	container_reference<self_type> (), it2_ () {}
2111	BOOST_UBLAS_INLINE
2112	iterator2 (self_type &m, const subiterator2_type &it2):
2113	container_reference<self_type> (m), it2_ (it2) {}
2114
2115	// Arithmetic
2116	BOOST_UBLAS_INLINE
2117	iterator2 &operator ++ () {
2118	++ it2_;
2119	return *this;
2120	}
2121	BOOST_UBLAS_INLINE
2122	iterator2 &operator -- () {
2123	-- it2_;
2124	return *this;
2125	}
2126	BOOST_UBLAS_INLINE
2127	iterator2 &operator += (difference_type n) {
2128	it2_ += n;
2129	return *this;
2130	}
2131	BOOST_UBLAS_INLINE
2132	iterator2 &operator -= (difference_type n) {
2133	it2_ -= n;
2134	return *this;
2135	}
2136	BOOST_UBLAS_INLINE
2137	difference_type operator - (const iterator2 &it) const {
2138	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2139	return it2_ - it.it2_;
2140	}
2141
2142	// Dereference
2143	BOOST_UBLAS_INLINE
2144	reference operator * () const {
2145	size_type i = index1 ();
2146	size_type j = index2 ();
2147	BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());
2148	BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());
2149	#ifdef BOOST_UBLAS_OWN_BANDED
2150	size_type k = (std::max) (i, j);
2151	size_type l = (*this) ().lower () + j - i;
2152	if (k < (std::max) ((*this) ().size1 (), (*this) ().size2 ()) &&
2153	l < (*this) ().lower () + 1 + (*this) ().upper ())
2154	return *it2_;
2155	#else
2156	size_type k = j;
2157	size_type l = (*this) ().upper () + i - j;
2158	if (k < (*this) ().size2 () &&
2159	l < (*this) ().lower () + 1 + (*this) ().upper ())
2160	return *it2_;
2161	#endif
2162	return (*this) () (i, j);
2163	}
2164	BOOST_UBLAS_INLINE
2165	reference operator [] (difference_type n) const {
2166	return *(*this + n);
2167	}
2168
2169	#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
2170	BOOST_UBLAS_INLINE
2171	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2172	typename self_type::
2173	#endif
2174	iterator1 begin () const {
2175	return (*this) ().find1 (1, 0, index2 ());
2176	}
2177	BOOST_UBLAS_INLINE
2178	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2179	typename self_type::
2180	#endif
2181	iterator1 end () const {
2182	return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
2183	}
2184	BOOST_UBLAS_INLINE
2185	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2186	typename self_type::
2187	#endif
2188	reverse_iterator1 rbegin () const {
2189	return reverse_iterator1 (end ());
2190	}
2191	BOOST_UBLAS_INLINE
2192	#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2193	typename self_type::
2194	#endif
2195	reverse_iterator1 rend () const {
2196	return reverse_iterator1 (begin ());
2197	}
2198	#endif
2199
2200	// Indices
2201	BOOST_UBLAS_INLINE
2202	size_type index1 () const {
2203	return it2_.index1 ();
2204	}
2205	BOOST_UBLAS_INLINE
2206	size_type index2 () const {
2207	return it2_.index2 ();
2208	}
2209
2210	// Assignment
2211	BOOST_UBLAS_INLINE
2212	iterator2 &operator = (const iterator2 &it) {
2213	container_reference<self_type>::assign (&it ());
2214	it2_ = it.it2_;
2215	return *this;
2216	}
2217
2218	// Comparison
2219	BOOST_UBLAS_INLINE
2220	bool operator == (const iterator2 &it) const {
2221	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2222	return it2_ == it.it2_;
2223	}
2224	BOOST_UBLAS_INLINE
2225	bool operator < (const iterator2 &it) const {
2226	BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
2227	return it2_ < it.it2_;
2228	}
2229
2230	private:
2231	subiterator2_type it2_;
2232
2233	friend class const_iterator2;
2234	};
2235	#endif
2236
2237	BOOST_UBLAS_INLINE
2238	iterator2 begin2 () {
2239	return find2 (0, 0, 0);
2240	}
2241	BOOST_UBLAS_INLINE
2242	iterator2 end2 () {
2243	return find2 (0, 0, size2 ());
2244	}
2245
2246	// Reverse iterators
2247
2248	BOOST_UBLAS_INLINE
2249	const_reverse_iterator1 rbegin1 () const {
2250	return const_reverse_iterator1 (end1 ());
2251	}
2252	BOOST_UBLAS_INLINE
2253	const_reverse_iterator1 crbegin1 () const {
2254	return rbegin1 ();
2255	}
2256	BOOST_UBLAS_INLINE
2257	const_reverse_iterator1 rend1 () const {
2258	return const_reverse_iterator1 (begin1 ());
2259	}
2260	BOOST_UBLAS_INLINE
2261	const_reverse_iterator1 crend1 () const {
2262	return rend1 ();
2263	}
2264
2265	BOOST_UBLAS_INLINE
2266	reverse_iterator1 rbegin1 () {
2267	return reverse_iterator1 (end1 ());
2268	}
2269	BOOST_UBLAS_INLINE
2270	reverse_iterator1 rend1 () {
2271	return reverse_iterator1 (begin1 ());
2272	}
2273
2274	BOOST_UBLAS_INLINE
2275	const_reverse_iterator2 rbegin2 () const {
2276	return const_reverse_iterator2 (end2 ());
2277	}
2278	BOOST_UBLAS_INLINE
2279	const_reverse_iterator2 crbegin2 () const {
2280	return rbegin2 ();
2281	}
2282	BOOST_UBLAS_INLINE
2283	const_reverse_iterator2 rend2 () const {
2284	return const_reverse_iterator2 (begin2 ());
2285	}
2286	BOOST_UBLAS_INLINE
2287	const_reverse_iterator2 crend2 () const {
2288	return rend2 ();
2289	}
2290
2291	BOOST_UBLAS_INLINE
2292	reverse_iterator2 rbegin2 () {
2293	return reverse_iterator2 (end2 ());
2294	}
2295	BOOST_UBLAS_INLINE
2296	reverse_iterator2 rend2 () {
2297	return reverse_iterator2 (begin2 ());
2298	}
2299
2300	private:
2301	matrix_closure_type data_;
2302	size_type lower_;
2303	size_type upper_;
2304	typedef const value_type const_value_type;
2305	static const_value_type zero_;
2306	};
2307
2308	// Specialization for temporary_traits
2309	template <class M>
2310	struct vector_temporary_traits< banded_adaptor<M> >
2311	: vector_temporary_traits< M > {} ;
2312	template <class M>
2313	struct vector_temporary_traits< const banded_adaptor<M> >
2314	: vector_temporary_traits< M > {} ;
2315
2316	template <class M>
2317	struct matrix_temporary_traits< banded_adaptor<M> >
2318	: matrix_temporary_traits< M > {} ;
2319	template <class M>
2320	struct matrix_temporary_traits< const banded_adaptor<M> >
2321	: matrix_temporary_traits< M > {} ;
2322
2323
2324	template<class M>
2325	typename banded_adaptor<M>::const_value_type banded_adaptor<M>::zero_ = value_type/*zero*/();
2326
2327	/** \brief A diagonal matrix adaptator: convert a any matrix into a diagonal matrix expression
2328	*
2329	* For a \f$(m\times m)\f$-dimensional matrix, the \c diagonal_adaptor will provide a diagonal matrix
2330	* with \f$0 \leq i < m\f$ and \f$0 \leq j < m\f$, if \f$i\neq j\f$ then \f$b_{i,j}=0\f$.
2331	*
2332	* Storage and location are based on those of the underlying matrix. This is important because
2333	* a \c diagonal_adaptor does not copy the matrix data to a new place. Therefore, modifying values
2334	* in a \c diagonal_adaptor matrix will also modify the underlying matrix too.
2335	*
2336	* \tparam M the type of matrix used to generate the diagonal matrix
2337	*/
2338
2339	template<class M>
2340	class diagonal_adaptor:
2341	public banded_adaptor<M> {
2342	public:
2343	typedef M matrix_type;
2344	typedef banded_adaptor<M> adaptor_type;
2345
2346	// Construction and destruction
2347	BOOST_UBLAS_INLINE
2348	diagonal_adaptor ():
2349	adaptor_type () {}
2350	BOOST_UBLAS_INLINE
2351	diagonal_adaptor (matrix_type &data):
2352	adaptor_type (data) {}
2353	BOOST_UBLAS_INLINE
2354	~diagonal_adaptor () {}
2355
2356	// Assignment
2357	BOOST_UBLAS_INLINE
2358	diagonal_adaptor &operator = (const diagonal_adaptor &m) {
2359	adaptor_type::operator = (m);
2360	return *this;
2361	}
2362	template<class AE>
2363	BOOST_UBLAS_INLINE
2364	diagonal_adaptor &operator = (const matrix_expression<AE> &ae) {
2365	adaptor_type::operator = (ae);
2366	return *this;
2367	}
2368	};
2369
2370	}}}
2371
2372	#endif
2373