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

Revision 857, 96.8 KB checked in by igarcia, 19 years ago (diff)
RevLine 
[857]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
26namespace 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.