Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

triangular.hpp @ 857

Revision 857, 96.8 KB checked in by igarcia, 18 years ago (diff)

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

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

Context Navigation

source: NonGTP/Boost/boost/numeric/ublas/triangular.hpp @ 857

Download in other formats: