Planet

navi

home

PPS

about

screenshots

download

development

forum

Context Navigation

source: downloads/boost_1_34_1/boost/numeric/ublas/triangular.hpp @ 33

Last change on this file since 33 was 29, checked in by landauf, 16 years ago
updated boost from 1_33_1 to 1_34_1
File size: 95.3 KB

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

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

Download in other formats: