Planet
navi homePPSaboutscreenshotsdownloaddevelopmentforum

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

Last change on this file since 47 was 29, checked in by landauf, 16 years ago

updated boost from 1_33_1 to 1_34_1

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