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