Planet

navi

home

PPS

about

screenshots

download

development

forum

Context Navigation

source: downloads/boost_1_33_1/boost/numeric/ublas/triangular.hpp @ 12

Last change on this file since 12 was 12, checked in by landauf, 17 years ago
added boost
File size: 94.3 KB

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

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

Download in other formats: