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