Symaptic: os/ossrv/ossrv_pub/boost_apis/boost/numeric/ublas/triangular.hpp@260cb5ec6c19 (annotated)

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

author	sl
	Tue, 10 Jun 2014 14:32:02 +0200
changeset 1	260cb5ec6c19
permissions	-rw-r--r--