Boost.uBlas 1.49
Linear Algebra in C++: matrices, vectors and numeric algorithms

triangular.hpp

Go to the documentation of this file.
00001 //
00002 //  Copyright (c) 2000-2002
00003 //  Joerg Walter, Mathias Koch
00004 //
00005 //  Distributed under the Boost Software License, Version 1.0. (See
00006 //  accompanying file LICENSE_1_0.txt or copy at
00007 //  http://www.boost.org/LICENSE_1_0.txt)
00008 //
00009 //  The authors gratefully acknowledge the support of
00010 //  GeNeSys mbH & Co. KG in producing this work.
00011 //
00012 
00013 #ifndef _BOOST_UBLAS_TRIANGULAR_
00014 #define _BOOST_UBLAS_TRIANGULAR_
00015 
00016 #include <boost/numeric/ublas/matrix.hpp>
00017 #include <boost/numeric/ublas/detail/temporary.hpp>
00018 #include <boost/type_traits/remove_const.hpp>
00019 
00020 // Iterators based on ideas of Jeremy Siek
00021 
00022 namespace boost { namespace numeric { namespace ublas {
00023 
00024     namespace detail {
00025         using namespace boost::numeric::ublas;
00026 
00027         // Matrix resizing algorithm
00028         template <class L, class T, class M>
00029         BOOST_UBLAS_INLINE
00030         void matrix_resize_preserve (M& m, M& temporary) {
00031             typedef L layout_type;
00032             typedef T triangular_type;
00033             typedef typename M::size_type size_type;
00034             const size_type msize1 (m.size1 ());        // original size
00035             const size_type msize2 (m.size2 ());
00036             const size_type size1 (temporary.size1 ());    // new size is specified by temporary
00037             const size_type size2 (temporary.size2 ());
00038             // Common elements to preserve
00039             const size_type size1_min = (std::min) (size1, msize1);
00040             const size_type size2_min = (std::min) (size2, msize2);
00041             // Order for major and minor sizes
00042             const size_type major_size = layout_type::size_M (size1_min, size2_min);
00043             const size_type minor_size = layout_type::size_m (size1_min, size2_min);
00044             // Indexing copy over major
00045             for (size_type major = 0; major != major_size; ++major) {
00046                 for (size_type minor = 0; minor != minor_size; ++minor) {
00047                         // find indexes - use invertability of element_ functions
00048                     const size_type i1 = layout_type::index_M(major, minor);
00049                     const size_type i2 = layout_type::index_m(major, minor);
00050                     if ( triangular_type::other(i1,i2) ) {
00051                         temporary.data () [triangular_type::element (layout_type (), i1, size1, i2, size2)] =
00052                             m.data() [triangular_type::element (layout_type (), i1, msize1, i2, msize2)];
00053                     }
00054                 }
00055             }
00056             m.assign_temporary (temporary);
00057         }
00058     }
00059 
00077     template<class T, class TRI, class L, class A>
00078     class triangular_matrix:
00079         public matrix_container<triangular_matrix<T, TRI, L, A> > {
00080 
00081         typedef T *pointer;
00082         typedef TRI triangular_type;
00083         typedef L layout_type;
00084         typedef triangular_matrix<T, TRI, L, A> self_type;
00085     public:
00086 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
00087         using matrix_container<self_type>::operator ();
00088 #endif
00089         typedef typename A::size_type size_type;
00090         typedef typename A::difference_type difference_type;
00091         typedef T value_type;
00092         typedef const T &const_reference;
00093         typedef T &reference;
00094         typedef A array_type;
00095 
00096         typedef const matrix_reference<const self_type> const_closure_type;
00097         typedef matrix_reference<self_type> closure_type;
00098         typedef vector<T, A> vector_temporary_type;
00099         typedef matrix<T, L, A> matrix_temporary_type;  // general sub-matrix
00100         typedef packed_tag storage_category;
00101         typedef typename L::orientation_category orientation_category;
00102 
00103         // Construction and destruction
00104         BOOST_UBLAS_INLINE
00105         triangular_matrix ():
00106             matrix_container<self_type> (),
00107             size1_ (0), size2_ (0), data_ (0) {}
00108         BOOST_UBLAS_INLINE
00109         triangular_matrix (size_type size1, size_type size2):
00110             matrix_container<self_type> (),
00111             size1_ (size1), size2_ (size2), data_ (triangular_type::packed_size (layout_type (), size1, size2)) {
00112         }
00113         BOOST_UBLAS_INLINE
00114         triangular_matrix (size_type size1, size_type size2, const array_type &data):
00115             matrix_container<self_type> (),
00116             size1_ (size1), size2_ (size2), data_ (data) {}
00117         BOOST_UBLAS_INLINE
00118         triangular_matrix (const triangular_matrix &m):
00119             matrix_container<self_type> (),
00120             size1_ (m.size1_), size2_ (m.size2_), data_ (m.data_) {}
00121         template<class AE>
00122         BOOST_UBLAS_INLINE
00123         triangular_matrix (const matrix_expression<AE> &ae):
00124             matrix_container<self_type> (),
00125             size1_ (ae ().size1 ()), size2_ (ae ().size2 ()),
00126             data_ (triangular_type::packed_size (layout_type (), size1_, size2_)) {
00127             matrix_assign<scalar_assign> (*this, ae);
00128         }
00129 
00130         // Accessors
00131         BOOST_UBLAS_INLINE
00132         size_type size1 () const {
00133             return size1_;
00134         }
00135         BOOST_UBLAS_INLINE
00136         size_type size2 () const {
00137             return size2_;
00138         }
00139 
00140         // Storage accessors
00141         BOOST_UBLAS_INLINE
00142         const array_type &data () const {
00143             return data_;
00144         }
00145         BOOST_UBLAS_INLINE
00146         array_type &data () {
00147             return data_;
00148         }
00149 
00150         // Resizing
00151         BOOST_UBLAS_INLINE
00152         void resize (size_type size1, size_type size2, bool preserve = true) {
00153             if (preserve) {
00154                 self_type temporary (size1, size2);
00155                 detail::matrix_resize_preserve<layout_type, triangular_type> (*this, temporary);
00156             }
00157             else {
00158                 data ().resize (triangular_type::packed_size (layout_type (), size1, size2));
00159                 size1_ = size1;
00160                 size2_ = size2;
00161             }
00162         }
00163         BOOST_UBLAS_INLINE
00164         void resize_packed_preserve (size_type size1, size_type size2) {
00165             size1_ = size1;
00166             size2_ = size2;
00167             data ().resize (triangular_type::packed_size (layout_type (), size1_, size2_), value_type ());
00168         }
00169 
00170         // Element access
00171         BOOST_UBLAS_INLINE
00172         const_reference operator () (size_type i, size_type j) const {
00173             BOOST_UBLAS_CHECK (i < size1_, bad_index ());
00174             BOOST_UBLAS_CHECK (j < size2_, bad_index ());
00175             if (triangular_type::other (i, j))
00176                 return data () [triangular_type::element (layout_type (), i, size1_, j, size2_)];
00177             else if (triangular_type::one (i, j))
00178                 return one_;
00179             else
00180                 return zero_;
00181         }
00182         BOOST_UBLAS_INLINE
00183         reference at_element (size_type i, size_type j) {
00184             BOOST_UBLAS_CHECK (i < size1_, bad_index ());
00185             BOOST_UBLAS_CHECK (j < size2_, bad_index ());
00186             return data () [triangular_type::element (layout_type (), i, size1_, j, size2_)];
00187         }
00188         BOOST_UBLAS_INLINE
00189         reference operator () (size_type i, size_type j) {
00190             BOOST_UBLAS_CHECK (i < size1_, bad_index ());
00191             BOOST_UBLAS_CHECK (j < size2_, bad_index ());
00192             if (!triangular_type::other (i, j)) {
00193                 bad_index ().raise ();
00194                 // NEVER reached
00195             }
00196             return data () [triangular_type::element (layout_type (), i, size1_, j, size2_)];
00197         }
00198         
00199         // Element assignment
00200         BOOST_UBLAS_INLINE
00201         reference insert_element (size_type i, size_type j, const_reference t) {
00202             return (operator () (i, j) = t);
00203         }
00204         BOOST_UBLAS_INLINE
00205         void erase_element (size_type i, size_type j) {
00206             operator () (i, j) = value_type/*zero*/();
00207         }
00208         
00209         // Zeroing
00210         BOOST_UBLAS_INLINE
00211         void clear () {
00212             // data ().clear ();
00213             std::fill (data ().begin (), data ().end (), value_type/*zero*/());
00214         }
00215 
00216         // Assignment
00217         BOOST_UBLAS_INLINE
00218         triangular_matrix &operator = (const triangular_matrix &m) {
00219             size1_ = m.size1_;
00220             size2_ = m.size2_;
00221             data () = m.data ();
00222             return *this;
00223         }
00224         BOOST_UBLAS_INLINE
00225         triangular_matrix &assign_temporary (triangular_matrix &m) {
00226             swap (m);
00227             return *this;
00228         }
00229         template<class AE>
00230         BOOST_UBLAS_INLINE
00231         triangular_matrix &operator = (const matrix_expression<AE> &ae) {
00232             self_type temporary (ae);
00233             return assign_temporary (temporary);
00234         }
00235         template<class AE>
00236         BOOST_UBLAS_INLINE
00237         triangular_matrix &assign (const matrix_expression<AE> &ae) {
00238             matrix_assign<scalar_assign> (*this, ae);
00239             return *this;
00240         }
00241         template<class AE>
00242         BOOST_UBLAS_INLINE
00243         triangular_matrix& operator += (const matrix_expression<AE> &ae) {
00244             self_type temporary (*this + ae);
00245             return assign_temporary (temporary);
00246         }
00247         template<class AE>
00248         BOOST_UBLAS_INLINE
00249         triangular_matrix &plus_assign (const matrix_expression<AE> &ae) {
00250             matrix_assign<scalar_plus_assign> (*this, ae);
00251             return *this;
00252         }
00253         template<class AE>
00254         BOOST_UBLAS_INLINE
00255         triangular_matrix& operator -= (const matrix_expression<AE> &ae) {
00256             self_type temporary (*this - ae);
00257             return assign_temporary (temporary);
00258         }
00259         template<class AE>
00260         BOOST_UBLAS_INLINE
00261         triangular_matrix &minus_assign (const matrix_expression<AE> &ae) {
00262             matrix_assign<scalar_minus_assign> (*this, ae);
00263             return *this;
00264         }
00265         template<class AT>
00266         BOOST_UBLAS_INLINE
00267         triangular_matrix& operator *= (const AT &at) {
00268             matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
00269             return *this;
00270         }
00271         template<class AT>
00272         BOOST_UBLAS_INLINE
00273         triangular_matrix& operator /= (const AT &at) {
00274             matrix_assign_scalar<scalar_divides_assign> (*this, at);
00275             return *this;
00276         }
00277 
00278         // Swapping
00279         BOOST_UBLAS_INLINE
00280         void swap (triangular_matrix &m) {
00281             if (this != &m) {
00282                 // BOOST_UBLAS_CHECK (size2_ == m.size2_, bad_size ());
00283                 std::swap (size1_, m.size1_);
00284                 std::swap (size2_, m.size2_);
00285                 data ().swap (m.data ());
00286             }
00287         }
00288         BOOST_UBLAS_INLINE
00289         friend void swap (triangular_matrix &m1, triangular_matrix &m2) {
00290             m1.swap (m2);
00291         }
00292 
00293         // Iterator types
00294 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
00295         typedef indexed_iterator1<self_type, packed_random_access_iterator_tag> iterator1;
00296         typedef indexed_iterator2<self_type, packed_random_access_iterator_tag> iterator2;
00297         typedef indexed_const_iterator1<self_type, packed_random_access_iterator_tag> const_iterator1;
00298         typedef indexed_const_iterator2<self_type, packed_random_access_iterator_tag> const_iterator2;
00299 #else
00300         class const_iterator1;
00301         class iterator1;
00302         class const_iterator2;
00303         class iterator2;
00304 #endif
00305         typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
00306         typedef reverse_iterator_base1<iterator1> reverse_iterator1;
00307         typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
00308         typedef reverse_iterator_base2<iterator2> reverse_iterator2;
00309 
00310         // Element lookup
00311         BOOST_UBLAS_INLINE
00312         const_iterator1 find1 (int rank, size_type i, size_type j) const {
00313             if (rank == 1)
00314                 i = triangular_type::restrict1 (i, j, size1_, size2_);
00315             if (rank == 0)
00316                 i = triangular_type::global_restrict1 (i, size1_, j, size2_);
00317             return const_iterator1 (*this, i, j);
00318         }
00319         BOOST_UBLAS_INLINE
00320         iterator1 find1 (int rank, size_type i, size_type j) {
00321             if (rank == 1)
00322                 i = triangular_type::mutable_restrict1 (i, j, size1_, size2_);
00323             if (rank == 0)
00324                 i = triangular_type::global_mutable_restrict1 (i, size1_, j, size2_);
00325             return iterator1 (*this, i, j);
00326         }
00327         BOOST_UBLAS_INLINE
00328         const_iterator2 find2 (int rank, size_type i, size_type j) const {
00329             if (rank == 1)
00330                 j = triangular_type::restrict2 (i, j, size1_, size2_);
00331             if (rank == 0)
00332                 j = triangular_type::global_restrict2 (i, size1_, j, size2_);
00333             return const_iterator2 (*this, i, j);
00334         }
00335         BOOST_UBLAS_INLINE
00336         iterator2 find2 (int rank, size_type i, size_type j) {
00337             if (rank == 1)
00338                 j = triangular_type::mutable_restrict2 (i, j, size1_, size2_);
00339             if (rank == 0)
00340                 j = triangular_type::global_mutable_restrict2 (i, size1_, j, size2_);
00341             return iterator2 (*this, i, j);
00342         }
00343 
00344         // Iterators simply are indices.
00345 
00346 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
00347         class const_iterator1:
00348             public container_const_reference<triangular_matrix>,
00349             public random_access_iterator_base<packed_random_access_iterator_tag,
00350                                                const_iterator1, value_type> {
00351         public:
00352             typedef typename triangular_matrix::value_type value_type;
00353             typedef typename triangular_matrix::difference_type difference_type;
00354             typedef typename triangular_matrix::const_reference reference;
00355             typedef const typename triangular_matrix::pointer pointer;
00356 
00357             typedef const_iterator2 dual_iterator_type;
00358             typedef const_reverse_iterator2 dual_reverse_iterator_type;
00359 
00360             // Construction and destruction
00361             BOOST_UBLAS_INLINE
00362             const_iterator1 ():
00363                 container_const_reference<self_type> (), it1_ (), it2_ () {}
00364             BOOST_UBLAS_INLINE
00365             const_iterator1 (const self_type &m, size_type it1, size_type it2):
00366                 container_const_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
00367             BOOST_UBLAS_INLINE
00368             const_iterator1 (const iterator1 &it):
00369                 container_const_reference<self_type> (it ()), it1_ (it.it1_), it2_ (it.it2_) {}
00370 
00371             // Arithmetic
00372             BOOST_UBLAS_INLINE
00373             const_iterator1 &operator ++ () {
00374                 ++ it1_;
00375                 return *this;
00376             }
00377             BOOST_UBLAS_INLINE
00378             const_iterator1 &operator -- () {
00379                 -- it1_;
00380                 return *this;
00381             }
00382             BOOST_UBLAS_INLINE
00383             const_iterator1 &operator += (difference_type n) {
00384                 it1_ += n;
00385                 return *this;
00386             }
00387             BOOST_UBLAS_INLINE
00388             const_iterator1 &operator -= (difference_type n) {
00389                 it1_ -= n;
00390                 return *this;
00391             }
00392             BOOST_UBLAS_INLINE
00393             difference_type operator - (const const_iterator1 &it) const {
00394                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
00395                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
00396                 return it1_ - it.it1_;
00397             }
00398 
00399             // Dereference
00400             BOOST_UBLAS_INLINE
00401             const_reference operator * () const {
00402                 return (*this) () (it1_, it2_);
00403             }
00404             BOOST_UBLAS_INLINE
00405             const_reference operator [] (difference_type n) const {
00406                 return *(*this + n);
00407             }
00408 
00409 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00410             BOOST_UBLAS_INLINE
00411 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
00412             typename self_type::
00413 #endif
00414             const_iterator2 begin () const {
00415                 return (*this) ().find2 (1, it1_, 0);
00416             }
00417             BOOST_UBLAS_INLINE
00418 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
00419             typename self_type::
00420 #endif
00421             const_iterator2 end () const {
00422                 return (*this) ().find2 (1, it1_, (*this) ().size2 ());
00423             }
00424             BOOST_UBLAS_INLINE
00425 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
00426             typename self_type::
00427 #endif
00428             const_reverse_iterator2 rbegin () const {
00429                 return const_reverse_iterator2 (end ());
00430             }
00431             BOOST_UBLAS_INLINE
00432 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
00433             typename self_type::
00434 #endif
00435             const_reverse_iterator2 rend () const {
00436                 return const_reverse_iterator2 (begin ());
00437             }
00438 #endif
00439 
00440             // Indices
00441             BOOST_UBLAS_INLINE
00442             size_type index1 () const {
00443                 return it1_;
00444             }
00445             BOOST_UBLAS_INLINE
00446             size_type index2 () const {
00447                 return it2_;
00448             }
00449 
00450             // Assignment
00451             BOOST_UBLAS_INLINE
00452             const_iterator1 &operator = (const const_iterator1 &it) {
00453                 container_const_reference<self_type>::assign (&it ());
00454                 it1_ = it.it1_;
00455                 it2_ = it.it2_;
00456                 return *this;
00457             }
00458 
00459             // Comparison
00460             BOOST_UBLAS_INLINE
00461             bool operator == (const const_iterator1 &it) const {
00462                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
00463                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
00464                 return it1_ == it.it1_;
00465             }
00466             BOOST_UBLAS_INLINE
00467             bool operator < (const const_iterator1 &it) const {
00468                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
00469                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
00470                 return it1_ < it.it1_;
00471             }
00472 
00473         private:
00474             size_type it1_;
00475             size_type it2_;
00476         };
00477 #endif
00478 
00479         BOOST_UBLAS_INLINE
00480         const_iterator1 begin1 () const {
00481             return find1 (0, 0, 0);
00482         }
00483         BOOST_UBLAS_INLINE
00484         const_iterator1 end1 () const {
00485             return find1 (0, size1_, 0);
00486         }
00487 
00488 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
00489         class iterator1:
00490             public container_reference<triangular_matrix>,
00491             public random_access_iterator_base<packed_random_access_iterator_tag,
00492                                                iterator1, value_type> {
00493         public:
00494             typedef typename triangular_matrix::value_type value_type;
00495             typedef typename triangular_matrix::difference_type difference_type;
00496             typedef typename triangular_matrix::reference reference;
00497             typedef typename triangular_matrix::pointer pointer;
00498 
00499             typedef iterator2 dual_iterator_type;
00500             typedef reverse_iterator2 dual_reverse_iterator_type;
00501 
00502             // Construction and destruction
00503             BOOST_UBLAS_INLINE
00504             iterator1 ():
00505                 container_reference<self_type> (), it1_ (), it2_ () {}
00506             BOOST_UBLAS_INLINE
00507             iterator1 (self_type &m, size_type it1, size_type it2):
00508                 container_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
00509 
00510             // Arithmetic
00511             BOOST_UBLAS_INLINE
00512             iterator1 &operator ++ () {
00513                 ++ it1_;
00514                 return *this;
00515             }
00516             BOOST_UBLAS_INLINE
00517             iterator1 &operator -- () {
00518                 -- it1_;
00519                 return *this;
00520             }
00521             BOOST_UBLAS_INLINE
00522             iterator1 &operator += (difference_type n) {
00523                 it1_ += n;
00524                 return *this;
00525             }
00526             BOOST_UBLAS_INLINE
00527             iterator1 &operator -= (difference_type n) {
00528                 it1_ -= n;
00529                 return *this;
00530             }
00531             BOOST_UBLAS_INLINE
00532             difference_type operator - (const iterator1 &it) const {
00533                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
00534                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
00535                 return it1_ - it.it1_;
00536             }
00537 
00538             // Dereference
00539             BOOST_UBLAS_INLINE
00540             reference operator * () const {
00541                 return (*this) () (it1_, it2_);
00542             }
00543             BOOST_UBLAS_INLINE
00544             reference operator [] (difference_type n) const {
00545                 return *(*this + n);
00546             }
00547 
00548 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00549             BOOST_UBLAS_INLINE
00550 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
00551             typename self_type::
00552 #endif
00553             iterator2 begin () const {
00554                 return (*this) ().find2 (1, it1_, 0);
00555             }
00556             BOOST_UBLAS_INLINE
00557 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
00558             typename self_type::
00559 #endif
00560             iterator2 end () const {
00561                 return (*this) ().find2 (1, it1_, (*this) ().size2 ());
00562             }
00563             BOOST_UBLAS_INLINE
00564 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
00565             typename self_type::
00566 #endif
00567             reverse_iterator2 rbegin () const {
00568                 return reverse_iterator2 (end ());
00569             }
00570             BOOST_UBLAS_INLINE
00571 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
00572             typename self_type::
00573 #endif
00574             reverse_iterator2 rend () const {
00575                 return reverse_iterator2 (begin ());
00576             }
00577 #endif
00578 
00579             // Indices
00580             BOOST_UBLAS_INLINE
00581             size_type index1 () const {
00582                 return it1_;
00583             }
00584             BOOST_UBLAS_INLINE
00585             size_type index2 () const {
00586                 return it2_;
00587             }
00588 
00589             // Assignment
00590             BOOST_UBLAS_INLINE
00591             iterator1 &operator = (const iterator1 &it) {
00592                 container_reference<self_type>::assign (&it ());
00593                 it1_ = it.it1_;
00594                 it2_ = it.it2_;
00595                 return *this;
00596             }
00597 
00598             // Comparison
00599             BOOST_UBLAS_INLINE
00600             bool operator == (const iterator1 &it) const {
00601                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
00602                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
00603                 return it1_ == it.it1_;
00604             }
00605             BOOST_UBLAS_INLINE
00606             bool operator < (const iterator1 &it) const {
00607                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
00608                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
00609                 return it1_ < it.it1_;
00610             }
00611 
00612         private:
00613             size_type it1_;
00614             size_type it2_;
00615 
00616             friend class const_iterator1;
00617         };
00618 #endif
00619 
00620         BOOST_UBLAS_INLINE
00621         iterator1 begin1 () {
00622             return find1 (0, 0, 0);
00623         }
00624         BOOST_UBLAS_INLINE
00625         iterator1 end1 () {
00626             return find1 (0, size1_, 0);
00627         }
00628 
00629 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
00630         class const_iterator2:
00631             public container_const_reference<triangular_matrix>,
00632             public random_access_iterator_base<packed_random_access_iterator_tag,
00633                                                const_iterator2, value_type> {
00634         public:
00635             typedef typename triangular_matrix::value_type value_type;
00636             typedef typename triangular_matrix::difference_type difference_type;
00637             typedef typename triangular_matrix::const_reference reference;
00638             typedef const typename triangular_matrix::pointer pointer;
00639 
00640             typedef const_iterator1 dual_iterator_type;
00641             typedef const_reverse_iterator1 dual_reverse_iterator_type;
00642 
00643             // Construction and destruction
00644             BOOST_UBLAS_INLINE
00645             const_iterator2 ():
00646                 container_const_reference<self_type> (), it1_ (), it2_ () {}
00647             BOOST_UBLAS_INLINE
00648             const_iterator2 (const self_type &m, size_type it1, size_type it2):
00649                 container_const_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
00650             BOOST_UBLAS_INLINE
00651             const_iterator2 (const iterator2 &it):
00652                 container_const_reference<self_type> (it ()), it1_ (it.it1_), it2_ (it.it2_) {}
00653 
00654             // Arithmetic
00655             BOOST_UBLAS_INLINE
00656             const_iterator2 &operator ++ () {
00657                 ++ it2_;
00658                 return *this;
00659             }
00660             BOOST_UBLAS_INLINE
00661             const_iterator2 &operator -- () {
00662                 -- it2_;
00663                 return *this;
00664             }
00665             BOOST_UBLAS_INLINE
00666             const_iterator2 &operator += (difference_type n) {
00667                 it2_ += n;
00668                 return *this;
00669             }
00670             BOOST_UBLAS_INLINE
00671             const_iterator2 &operator -= (difference_type n) {
00672                 it2_ -= n;
00673                 return *this;
00674             }
00675             BOOST_UBLAS_INLINE
00676             difference_type operator - (const const_iterator2 &it) const {
00677                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
00678                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
00679                 return it2_ - it.it2_;
00680             }
00681 
00682             // Dereference
00683             BOOST_UBLAS_INLINE
00684             const_reference operator * () const {
00685                 return (*this) () (it1_, it2_);
00686             }
00687             BOOST_UBLAS_INLINE
00688             const_reference operator [] (difference_type n) const {
00689                 return *(*this + n);
00690             }
00691 
00692 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00693             BOOST_UBLAS_INLINE
00694 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
00695             typename self_type::
00696 #endif
00697             const_iterator1 begin () const {
00698                 return (*this) ().find1 (1, 0, it2_);
00699             }
00700             BOOST_UBLAS_INLINE
00701 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
00702             typename self_type::
00703 #endif
00704             const_iterator1 end () const {
00705                 return (*this) ().find1 (1, (*this) ().size1 (), it2_);
00706             }
00707             BOOST_UBLAS_INLINE
00708 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
00709             typename self_type::
00710 #endif
00711             const_reverse_iterator1 rbegin () const {
00712                 return const_reverse_iterator1 (end ());
00713             }
00714             BOOST_UBLAS_INLINE
00715 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
00716             typename self_type::
00717 #endif
00718             const_reverse_iterator1 rend () const {
00719                 return const_reverse_iterator1 (begin ());
00720             }
00721 #endif
00722 
00723             // Indices
00724             BOOST_UBLAS_INLINE
00725             size_type index1 () const {
00726                 return it1_;
00727             }
00728             BOOST_UBLAS_INLINE
00729             size_type index2 () const {
00730                 return it2_;
00731             }
00732 
00733             // Assignment
00734             BOOST_UBLAS_INLINE
00735             const_iterator2 &operator = (const const_iterator2 &it) {
00736                 container_const_reference<self_type>::assign (&it ());
00737                 it1_ = it.it1_;
00738                 it2_ = it.it2_;
00739                 return *this;
00740             }
00741 
00742             // Comparison
00743             BOOST_UBLAS_INLINE
00744             bool operator == (const const_iterator2 &it) const {
00745                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
00746                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
00747                 return it2_ == it.it2_;
00748             }
00749             BOOST_UBLAS_INLINE
00750             bool operator < (const const_iterator2 &it) const {
00751                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
00752                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
00753                 return it2_ < it.it2_;
00754             }
00755 
00756         private:
00757             size_type it1_;
00758             size_type it2_;
00759         };
00760 #endif
00761 
00762         BOOST_UBLAS_INLINE
00763         const_iterator2 begin2 () const {
00764             return find2 (0, 0, 0);
00765         }
00766         BOOST_UBLAS_INLINE
00767         const_iterator2 end2 () const {
00768             return find2 (0, 0, size2_);
00769         }
00770 
00771 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
00772         class iterator2:
00773             public container_reference<triangular_matrix>,
00774             public random_access_iterator_base<packed_random_access_iterator_tag,
00775                                                iterator2, value_type> {
00776         public:
00777             typedef typename triangular_matrix::value_type value_type;
00778             typedef typename triangular_matrix::difference_type difference_type;
00779             typedef typename triangular_matrix::reference reference;
00780             typedef typename triangular_matrix::pointer pointer;
00781 
00782             typedef iterator1 dual_iterator_type;
00783             typedef reverse_iterator1 dual_reverse_iterator_type;
00784 
00785             // Construction and destruction
00786             BOOST_UBLAS_INLINE
00787             iterator2 ():
00788                 container_reference<self_type> (), it1_ (), it2_ () {}
00789             BOOST_UBLAS_INLINE
00790             iterator2 (self_type &m, size_type it1, size_type it2):
00791                 container_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
00792 
00793             // Arithmetic
00794             BOOST_UBLAS_INLINE
00795             iterator2 &operator ++ () {
00796                 ++ it2_;
00797                 return *this;
00798             }
00799             BOOST_UBLAS_INLINE
00800             iterator2 &operator -- () {
00801                 -- it2_;
00802                 return *this;
00803             }
00804             BOOST_UBLAS_INLINE
00805             iterator2 &operator += (difference_type n) {
00806                 it2_ += n;
00807                 return *this;
00808             }
00809             BOOST_UBLAS_INLINE
00810             iterator2 &operator -= (difference_type n) {
00811                 it2_ -= n;
00812                 return *this;
00813             }
00814             BOOST_UBLAS_INLINE
00815             difference_type operator - (const iterator2 &it) const {
00816                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
00817                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
00818                 return it2_ - it.it2_;
00819             }
00820 
00821             // Dereference
00822             BOOST_UBLAS_INLINE
00823             reference operator * () const {
00824                 return (*this) () (it1_, it2_);
00825             }
00826             BOOST_UBLAS_INLINE
00827             reference operator [] (difference_type n) const {
00828                 return *(*this + n);
00829             }
00830 
00831 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00832             BOOST_UBLAS_INLINE
00833 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
00834             typename self_type::
00835 #endif
00836             iterator1 begin () const {
00837                 return (*this) ().find1 (1, 0, it2_);
00838             }
00839             BOOST_UBLAS_INLINE
00840 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
00841             typename self_type::
00842 #endif
00843             iterator1 end () const {
00844                 return (*this) ().find1 (1, (*this) ().size1 (), it2_);
00845             }
00846             BOOST_UBLAS_INLINE
00847 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
00848             typename self_type::
00849 #endif
00850             reverse_iterator1 rbegin () const {
00851                 return reverse_iterator1 (end ());
00852             }
00853             BOOST_UBLAS_INLINE
00854 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
00855             typename self_type::
00856 #endif
00857             reverse_iterator1 rend () const {
00858                 return reverse_iterator1 (begin ());
00859             }
00860 #endif
00861 
00862             // Indices
00863             BOOST_UBLAS_INLINE
00864             size_type index1 () const {
00865                 return it1_;
00866             }
00867             BOOST_UBLAS_INLINE
00868             size_type index2 () const {
00869                 return it2_;
00870             }
00871 
00872             // Assignment
00873             BOOST_UBLAS_INLINE
00874             iterator2 &operator = (const iterator2 &it) {
00875                 container_reference<self_type>::assign (&it ());
00876                 it1_ = it.it1_;
00877                 it2_ = it.it2_;
00878                 return *this;
00879             }
00880 
00881             // Comparison
00882             BOOST_UBLAS_INLINE
00883             bool operator == (const iterator2 &it) const {
00884                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
00885                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
00886                 return it2_ == it.it2_;
00887             }
00888             BOOST_UBLAS_INLINE
00889             bool operator < (const iterator2 &it) const {
00890                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
00891                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
00892                 return it2_ < it.it2_;
00893             }
00894 
00895         private:
00896             size_type it1_;
00897             size_type it2_;
00898 
00899             friend class const_iterator2;
00900         };
00901 #endif
00902 
00903         BOOST_UBLAS_INLINE
00904         iterator2 begin2 () {
00905             return find2 (0, 0, 0);
00906         }
00907         BOOST_UBLAS_INLINE
00908         iterator2 end2 () {
00909             return find2 (0, 0, size2_);
00910         }
00911 
00912         // Reverse iterators
00913 
00914         BOOST_UBLAS_INLINE
00915         const_reverse_iterator1 rbegin1 () const {
00916             return const_reverse_iterator1 (end1 ());
00917         }
00918         BOOST_UBLAS_INLINE
00919         const_reverse_iterator1 rend1 () const {
00920             return const_reverse_iterator1 (begin1 ());
00921         }
00922 
00923         BOOST_UBLAS_INLINE
00924         reverse_iterator1 rbegin1 () {
00925             return reverse_iterator1 (end1 ());
00926         }
00927         BOOST_UBLAS_INLINE
00928         reverse_iterator1 rend1 () {
00929             return reverse_iterator1 (begin1 ());
00930         }
00931 
00932         BOOST_UBLAS_INLINE
00933         const_reverse_iterator2 rbegin2 () const {
00934             return const_reverse_iterator2 (end2 ());
00935         }
00936         BOOST_UBLAS_INLINE
00937         const_reverse_iterator2 rend2 () const {
00938             return const_reverse_iterator2 (begin2 ());
00939         }
00940 
00941         BOOST_UBLAS_INLINE
00942         reverse_iterator2 rbegin2 () {
00943             return reverse_iterator2 (end2 ());
00944         }
00945         BOOST_UBLAS_INLINE
00946         reverse_iterator2 rend2 () {
00947             return reverse_iterator2 (begin2 ());
00948         }
00949 
00950     private:
00951         size_type size1_;
00952         size_type size2_;
00953         array_type data_;
00954         static const value_type zero_;
00955         static const value_type one_;
00956     };
00957 
00958     template<class T, class TRI, class L, class A>
00959     const typename triangular_matrix<T, TRI, L, A>::value_type triangular_matrix<T, TRI, L, A>::zero_ = value_type/*zero*/();
00960     template<class T, class TRI, class L, class A>
00961     const typename triangular_matrix<T, TRI, L, A>::value_type triangular_matrix<T, TRI, L, A>::one_ (1);
00962 
00963 
00964     // Triangular matrix adaptor class
00965     template<class M, class TRI>
00966     class triangular_adaptor:
00967         public matrix_expression<triangular_adaptor<M, TRI> > {
00968 
00969         typedef triangular_adaptor<M, TRI> self_type;
00970 
00971     public:
00972 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
00973         using matrix_expression<self_type>::operator ();
00974 #endif
00975         typedef const M const_matrix_type;
00976         typedef M matrix_type;
00977         typedef TRI triangular_type;
00978         typedef typename M::size_type size_type;
00979         typedef typename M::difference_type difference_type;
00980         typedef typename M::value_type value_type;
00981         typedef typename M::const_reference const_reference;
00982         typedef typename boost::mpl::if_<boost::is_const<M>,
00983                                           typename M::const_reference,
00984                                           typename M::reference>::type reference;
00985         typedef typename boost::mpl::if_<boost::is_const<M>,
00986                                           typename M::const_closure_type,
00987                                           typename M::closure_type>::type matrix_closure_type;
00988         typedef const self_type const_closure_type;
00989         typedef self_type closure_type;
00990         // Replaced by _temporary_traits to avoid type requirements on M
00991         //typedef typename M::vector_temporary_type vector_temporary_type;
00992         //typedef typename M::matrix_temporary_type matrix_temporary_type;
00993         typedef typename storage_restrict_traits<typename M::storage_category,
00994                                                  packed_proxy_tag>::storage_category storage_category;
00995         typedef typename M::orientation_category orientation_category;
00996 
00997         // Construction and destruction
00998         BOOST_UBLAS_INLINE
00999         triangular_adaptor (matrix_type &data):
01000             matrix_expression<self_type> (),
01001             data_ (data) {}
01002         BOOST_UBLAS_INLINE
01003         triangular_adaptor (const triangular_adaptor &m):
01004             matrix_expression<self_type> (),
01005             data_ (m.data_) {}
01006 
01007         // Accessors
01008         BOOST_UBLAS_INLINE
01009         size_type size1 () const {
01010             return data_.size1 ();
01011         }
01012         BOOST_UBLAS_INLINE
01013         size_type size2 () const {
01014             return data_.size2 ();
01015         }
01016 
01017         // Storage accessors
01018         BOOST_UBLAS_INLINE
01019         const matrix_closure_type &data () const {
01020             return data_;
01021         }
01022         BOOST_UBLAS_INLINE
01023         matrix_closure_type &data () {
01024             return data_;
01025         }
01026 
01027         // Element access
01028 #ifndef BOOST_UBLAS_PROXY_CONST_MEMBER
01029         BOOST_UBLAS_INLINE
01030         const_reference operator () (size_type i, size_type j) const {
01031             BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
01032             BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
01033             if (triangular_type::other (i, j))
01034                 return data () (i, j);
01035             else if (triangular_type::one (i, j))
01036                 return one_;
01037             else
01038                 return zero_;
01039         }
01040         BOOST_UBLAS_INLINE
01041         reference operator () (size_type i, size_type j) {
01042             BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
01043             BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
01044             if (!triangular_type::other (i, j)) {
01045                 bad_index ().raise ();
01046                 // NEVER reached
01047             }
01048             return data () (i, j);
01049         }
01050 #else
01051         BOOST_UBLAS_INLINE
01052         reference operator () (size_type i, size_type j) const {
01053             BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
01054             BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
01055             if (!triangular_type::other (i, j)) {
01056                 bad_index ().raise ();
01057                 // NEVER reached
01058             }
01059             return data () (i, j);
01060         }
01061 #endif
01062 
01063         // Assignment
01064         BOOST_UBLAS_INLINE
01065         triangular_adaptor &operator = (const triangular_adaptor &m) {
01066             matrix_assign<scalar_assign> (*this, m);
01067             return *this;
01068         }
01069         BOOST_UBLAS_INLINE
01070         triangular_adaptor &assign_temporary (triangular_adaptor &m) {
01071             *this = m;
01072             return *this;
01073         }
01074         template<class AE>
01075         BOOST_UBLAS_INLINE
01076         triangular_adaptor &operator = (const matrix_expression<AE> &ae) {
01077             matrix_assign<scalar_assign> (*this, matrix<value_type> (ae));
01078             return *this;
01079         }
01080         template<class AE>
01081         BOOST_UBLAS_INLINE
01082         triangular_adaptor &assign (const matrix_expression<AE> &ae) {
01083             matrix_assign<scalar_assign> (*this, ae);
01084             return *this;
01085         }
01086         template<class AE>
01087         BOOST_UBLAS_INLINE
01088         triangular_adaptor& operator += (const matrix_expression<AE> &ae) {
01089             matrix_assign<scalar_assign> (*this, matrix<value_type> (*this + ae));
01090             return *this;
01091         }
01092         template<class AE>
01093         BOOST_UBLAS_INLINE
01094         triangular_adaptor &plus_assign (const matrix_expression<AE> &ae) {
01095             matrix_assign<scalar_plus_assign> (*this, ae);
01096             return *this;
01097         }
01098         template<class AE>
01099         BOOST_UBLAS_INLINE
01100         triangular_adaptor& operator -= (const matrix_expression<AE> &ae) {
01101             matrix_assign<scalar_assign> (*this, matrix<value_type> (*this - ae));
01102             return *this;
01103         }
01104         template<class AE>
01105         BOOST_UBLAS_INLINE
01106         triangular_adaptor &minus_assign (const matrix_expression<AE> &ae) {
01107             matrix_assign<scalar_minus_assign> (*this, ae);
01108             return *this;
01109         }
01110         template<class AT>
01111         BOOST_UBLAS_INLINE
01112         triangular_adaptor& operator *= (const AT &at) {
01113             matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
01114             return *this;
01115         }
01116         template<class AT>
01117         BOOST_UBLAS_INLINE
01118         triangular_adaptor& operator /= (const AT &at) {
01119             matrix_assign_scalar<scalar_divides_assign> (*this, at);
01120             return *this;
01121         }
01122 
01123         // Closure comparison
01124         BOOST_UBLAS_INLINE
01125         bool same_closure (const triangular_adaptor &ta) const {
01126             return (*this).data ().same_closure (ta.data ());
01127         }
01128 
01129         // Swapping
01130         BOOST_UBLAS_INLINE
01131         void swap (triangular_adaptor &m) {
01132             if (this != &m)
01133                 matrix_swap<scalar_swap> (*this, m);
01134         }
01135         BOOST_UBLAS_INLINE
01136         friend void swap (triangular_adaptor &m1, triangular_adaptor &m2) {
01137             m1.swap (m2);
01138         }
01139 
01140         // Iterator types
01141    private:
01142         typedef typename M::const_iterator1 const_subiterator1_type;
01143         typedef typename boost::mpl::if_<boost::is_const<M>,
01144                                           typename M::const_iterator1,
01145                                           typename M::iterator1>::type subiterator1_type;
01146         typedef typename M::const_iterator2 const_subiterator2_type;
01147         typedef typename boost::mpl::if_<boost::is_const<M>,
01148                                           typename M::const_iterator2,
01149                                           typename M::iterator2>::type subiterator2_type;
01150 
01151     public:
01152 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
01153         typedef indexed_iterator1<self_type, packed_random_access_iterator_tag> iterator1;
01154         typedef indexed_iterator2<self_type, packed_random_access_iterator_tag> iterator2;
01155         typedef indexed_const_iterator1<self_type, packed_random_access_iterator_tag> const_iterator1;
01156         typedef indexed_const_iterator2<self_type, packed_random_access_iterator_tag> const_iterator2;
01157 #else
01158         class const_iterator1;
01159         class iterator1;
01160         class const_iterator2;
01161         class iterator2;
01162 #endif
01163         typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
01164         typedef reverse_iterator_base1<iterator1> reverse_iterator1;
01165         typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
01166         typedef reverse_iterator_base2<iterator2> reverse_iterator2;
01167 
01168         // Element lookup
01169         BOOST_UBLAS_INLINE
01170         const_iterator1 find1 (int rank, size_type i, size_type j) const {
01171             if (rank == 1)
01172                 i = triangular_type::restrict1 (i, j, size1(), size2());
01173             if (rank == 0)
01174                 i = triangular_type::global_restrict1 (i, size1(), j, size2());
01175             return const_iterator1 (*this, data ().find1 (rank, i, j));
01176         }
01177         BOOST_UBLAS_INLINE
01178         iterator1 find1 (int rank, size_type i, size_type j) {
01179             if (rank == 1)
01180                 i = triangular_type::mutable_restrict1 (i, j, size1(), size2());
01181             if (rank == 0)
01182                 i = triangular_type::global_mutable_restrict1 (i, size1(), j, size2());
01183             return iterator1 (*this, data ().find1 (rank, i, j));
01184         }
01185         BOOST_UBLAS_INLINE
01186         const_iterator2 find2 (int rank, size_type i, size_type j) const {
01187             if (rank == 1)
01188                 j = triangular_type::restrict2 (i, j, size1(), size2());
01189             if (rank == 0)
01190                 j = triangular_type::global_restrict2 (i, size1(), j, size2());
01191             return const_iterator2 (*this, data ().find2 (rank, i, j));
01192         }
01193         BOOST_UBLAS_INLINE
01194         iterator2 find2 (int rank, size_type i, size_type j) {
01195             if (rank == 1)
01196                 j = triangular_type::mutable_restrict2 (i, j, size1(), size2());
01197             if (rank == 0)
01198                 j = triangular_type::global_mutable_restrict2 (i, size1(), j, size2());
01199             return iterator2 (*this, data ().find2 (rank, i, j));
01200         }
01201 
01202         // Iterators simply are indices.
01203 
01204 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
01205         class const_iterator1:
01206             public container_const_reference<triangular_adaptor>,
01207             public random_access_iterator_base<typename iterator_restrict_traits<
01208                                                    typename const_subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,
01209                                                const_iterator1, value_type> {
01210         public:
01211             typedef typename const_subiterator1_type::value_type value_type;
01212             typedef typename const_subiterator1_type::difference_type difference_type;
01213             typedef typename const_subiterator1_type::reference reference;
01214             typedef typename const_subiterator1_type::pointer pointer;
01215 
01216             typedef const_iterator2 dual_iterator_type;
01217             typedef const_reverse_iterator2 dual_reverse_iterator_type;
01218 
01219             // Construction and destruction
01220             BOOST_UBLAS_INLINE
01221             const_iterator1 ():
01222                 container_const_reference<self_type> (), it1_ () {}
01223             BOOST_UBLAS_INLINE
01224             const_iterator1 (const self_type &m, const const_subiterator1_type &it1):
01225                 container_const_reference<self_type> (m), it1_ (it1) {}
01226             BOOST_UBLAS_INLINE
01227             const_iterator1 (const iterator1 &it):
01228                 container_const_reference<self_type> (it ()), it1_ (it.it1_) {}
01229 
01230             // Arithmetic
01231             BOOST_UBLAS_INLINE
01232             const_iterator1 &operator ++ () {
01233                 ++ it1_;
01234                 return *this;
01235             }
01236             BOOST_UBLAS_INLINE
01237             const_iterator1 &operator -- () {
01238                 -- it1_;
01239                 return *this;
01240             }
01241             BOOST_UBLAS_INLINE
01242             const_iterator1 &operator += (difference_type n) {
01243                 it1_ += n;
01244                 return *this;
01245             }
01246             BOOST_UBLAS_INLINE
01247             const_iterator1 &operator -= (difference_type n) {
01248                 it1_ -= n;
01249                 return *this;
01250             }
01251             BOOST_UBLAS_INLINE
01252             difference_type operator - (const const_iterator1 &it) const {
01253                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
01254                 return it1_ - it.it1_;
01255             }
01256 
01257             // Dereference
01258             BOOST_UBLAS_INLINE
01259             const_reference operator * () const {
01260                 size_type i = index1 ();
01261                 size_type j = index2 ();
01262                 BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());
01263                 BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());
01264                 if (triangular_type::other (i, j))
01265                     return *it1_;
01266                 else
01267                     return (*this) () (i, j);
01268             }
01269             BOOST_UBLAS_INLINE
01270             const_reference operator [] (difference_type n) const {
01271                 return *(*this + n);
01272             }
01273 
01274 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
01275             BOOST_UBLAS_INLINE
01276 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
01277             typename self_type::
01278 #endif
01279             const_iterator2 begin () const {
01280                 return (*this) ().find2 (1, index1 (), 0);
01281             }
01282             BOOST_UBLAS_INLINE
01283 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
01284             typename self_type::
01285 #endif
01286             const_iterator2 end () const {
01287                 return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
01288             }
01289             BOOST_UBLAS_INLINE
01290 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
01291             typename self_type::
01292 #endif
01293             const_reverse_iterator2 rbegin () const {
01294                 return const_reverse_iterator2 (end ());
01295             }
01296             BOOST_UBLAS_INLINE
01297 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
01298             typename self_type::
01299 #endif
01300             const_reverse_iterator2 rend () const {
01301                 return const_reverse_iterator2 (begin ());
01302             }
01303 #endif
01304 
01305             // Indices
01306             BOOST_UBLAS_INLINE
01307             size_type index1 () const {
01308                 return it1_.index1 ();
01309             }
01310             BOOST_UBLAS_INLINE
01311             size_type index2 () const {
01312                 return it1_.index2 ();
01313             }
01314 
01315             // Assignment
01316             BOOST_UBLAS_INLINE
01317             const_iterator1 &operator = (const const_iterator1 &it) {
01318                 container_const_reference<self_type>::assign (&it ());
01319                 it1_ = it.it1_;
01320                 return *this;
01321             }
01322 
01323             // Comparison
01324             BOOST_UBLAS_INLINE
01325             bool operator == (const const_iterator1 &it) const {
01326                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
01327                 return it1_ == it.it1_;
01328             }
01329             BOOST_UBLAS_INLINE
01330             bool operator < (const const_iterator1 &it) const {
01331                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
01332                 return it1_ < it.it1_;
01333             }
01334 
01335         private:
01336             const_subiterator1_type it1_;
01337         };
01338 #endif
01339 
01340         BOOST_UBLAS_INLINE
01341         const_iterator1 begin1 () const {
01342             return find1 (0, 0, 0);
01343         }
01344         BOOST_UBLAS_INLINE
01345         const_iterator1 end1 () const {
01346             return find1 (0, size1 (), 0);
01347         }
01348 
01349 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
01350         class iterator1:
01351             public container_reference<triangular_adaptor>,
01352             public random_access_iterator_base<typename iterator_restrict_traits<
01353                                                    typename subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,
01354                                                iterator1, value_type> {
01355         public:
01356             typedef typename subiterator1_type::value_type value_type;
01357             typedef typename subiterator1_type::difference_type difference_type;
01358             typedef typename subiterator1_type::reference reference;
01359             typedef typename subiterator1_type::pointer pointer;
01360 
01361             typedef iterator2 dual_iterator_type;
01362             typedef reverse_iterator2 dual_reverse_iterator_type;
01363 
01364             // Construction and destruction
01365             BOOST_UBLAS_INLINE
01366             iterator1 ():
01367                 container_reference<self_type> (), it1_ () {}
01368             BOOST_UBLAS_INLINE
01369             iterator1 (self_type &m, const subiterator1_type &it1):
01370                 container_reference<self_type> (m), it1_ (it1) {}
01371 
01372             // Arithmetic
01373             BOOST_UBLAS_INLINE
01374             iterator1 &operator ++ () {
01375                 ++ it1_;
01376                 return *this;
01377             }
01378             BOOST_UBLAS_INLINE
01379             iterator1 &operator -- () {
01380                 -- it1_;
01381                 return *this;
01382             }
01383             BOOST_UBLAS_INLINE
01384             iterator1 &operator += (difference_type n) {
01385                 it1_ += n;
01386                 return *this;
01387             }
01388             BOOST_UBLAS_INLINE
01389             iterator1 &operator -= (difference_type n) {
01390                 it1_ -= n;
01391                 return *this;
01392             }
01393             BOOST_UBLAS_INLINE
01394             difference_type operator - (const iterator1 &it) const {
01395                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
01396                 return it1_ - it.it1_;
01397             }
01398 
01399             // Dereference
01400             BOOST_UBLAS_INLINE
01401             reference operator * () const {
01402                 size_type i = index1 ();
01403                 size_type j = index2 ();
01404                 BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());
01405                 BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());
01406                 if (triangular_type::other (i, j))
01407                     return *it1_;
01408                 else
01409                     return (*this) () (i, j);
01410             }
01411             BOOST_UBLAS_INLINE
01412             reference operator [] (difference_type n) const {
01413                 return *(*this + n);
01414             }
01415 
01416 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
01417             BOOST_UBLAS_INLINE
01418 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
01419             typename self_type::
01420 #endif
01421             iterator2 begin () const {
01422                 return (*this) ().find2 (1, index1 (), 0);
01423             }
01424             BOOST_UBLAS_INLINE
01425 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
01426             typename self_type::
01427 #endif
01428             iterator2 end () const {
01429                 return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
01430             }
01431             BOOST_UBLAS_INLINE
01432 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
01433             typename self_type::
01434 #endif
01435             reverse_iterator2 rbegin () const {
01436                 return reverse_iterator2 (end ());
01437             }
01438             BOOST_UBLAS_INLINE
01439 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
01440             typename self_type::
01441 #endif
01442             reverse_iterator2 rend () const {
01443                 return reverse_iterator2 (begin ());
01444             }
01445 #endif
01446 
01447             // Indices
01448             BOOST_UBLAS_INLINE
01449             size_type index1 () const {
01450                 return it1_.index1 ();
01451             }
01452             BOOST_UBLAS_INLINE
01453             size_type index2 () const {
01454                 return it1_.index2 ();
01455             }
01456 
01457             // Assignment
01458             BOOST_UBLAS_INLINE
01459             iterator1 &operator = (const iterator1 &it) {
01460                 container_reference<self_type>::assign (&it ());
01461                 it1_ = it.it1_;
01462                 return *this;
01463             }
01464 
01465             // Comparison
01466             BOOST_UBLAS_INLINE
01467             bool operator == (const iterator1 &it) const {
01468                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
01469                 return it1_ == it.it1_;
01470             }
01471             BOOST_UBLAS_INLINE
01472             bool operator < (const iterator1 &it) const {
01473                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
01474                 return it1_ < it.it1_;
01475             }
01476 
01477         private:
01478             subiterator1_type it1_;
01479 
01480             friend class const_iterator1;
01481         };
01482 #endif
01483 
01484         BOOST_UBLAS_INLINE
01485         iterator1 begin1 () {
01486             return find1 (0, 0, 0);
01487         }
01488         BOOST_UBLAS_INLINE
01489         iterator1 end1 () {
01490             return find1 (0, size1 (), 0);
01491         }
01492 
01493 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
01494         class const_iterator2:
01495             public container_const_reference<triangular_adaptor>,
01496             public random_access_iterator_base<typename iterator_restrict_traits<
01497                                                    typename const_subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,
01498                                                const_iterator2, value_type> {
01499         public:
01500             typedef typename const_subiterator2_type::value_type value_type;
01501             typedef typename const_subiterator2_type::difference_type difference_type;
01502             typedef typename const_subiterator2_type::reference reference;
01503             typedef typename const_subiterator2_type::pointer pointer;
01504 
01505             typedef const_iterator1 dual_iterator_type;
01506             typedef const_reverse_iterator1 dual_reverse_iterator_type;
01507 
01508             // Construction and destruction
01509             BOOST_UBLAS_INLINE
01510             const_iterator2 ():
01511                 container_const_reference<self_type> (), it2_ () {}
01512             BOOST_UBLAS_INLINE
01513             const_iterator2 (const self_type &m, const const_subiterator2_type &it2):
01514                 container_const_reference<self_type> (m), it2_ (it2) {}
01515             BOOST_UBLAS_INLINE
01516             const_iterator2 (const iterator2 &it):
01517                 container_const_reference<self_type> (it ()), it2_ (it.it2_) {}
01518 
01519             // Arithmetic
01520             BOOST_UBLAS_INLINE
01521             const_iterator2 &operator ++ () {
01522                 ++ it2_;
01523                 return *this;
01524             }
01525             BOOST_UBLAS_INLINE
01526             const_iterator2 &operator -- () {
01527                 -- it2_;
01528                 return *this;
01529             }
01530             BOOST_UBLAS_INLINE
01531             const_iterator2 &operator += (difference_type n) {
01532                 it2_ += n;
01533                 return *this;
01534             }
01535             BOOST_UBLAS_INLINE
01536             const_iterator2 &operator -= (difference_type n) {
01537                 it2_ -= n;
01538                 return *this;
01539             }
01540             BOOST_UBLAS_INLINE
01541             difference_type operator - (const const_iterator2 &it) const {
01542                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
01543                 return it2_ - it.it2_;
01544             }
01545 
01546             // Dereference
01547             BOOST_UBLAS_INLINE
01548             const_reference operator * () const {
01549                 size_type i = index1 ();
01550                 size_type j = index2 ();
01551                 BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());
01552                 BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());
01553                 if (triangular_type::other (i, j))
01554                     return *it2_;
01555                 else
01556                     return (*this) () (i, j);
01557             }
01558             BOOST_UBLAS_INLINE
01559             const_reference operator [] (difference_type n) const {
01560                 return *(*this + n);
01561             }
01562 
01563 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
01564             BOOST_UBLAS_INLINE
01565 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
01566             typename self_type::
01567 #endif
01568             const_iterator1 begin () const {
01569                 return (*this) ().find1 (1, 0, index2 ());
01570             }
01571             BOOST_UBLAS_INLINE
01572 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
01573             typename self_type::
01574 #endif
01575             const_iterator1 end () const {
01576                 return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
01577             }
01578             BOOST_UBLAS_INLINE
01579 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
01580             typename self_type::
01581 #endif
01582             const_reverse_iterator1 rbegin () const {
01583                 return const_reverse_iterator1 (end ());
01584             }
01585             BOOST_UBLAS_INLINE
01586 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
01587             typename self_type::
01588 #endif
01589             const_reverse_iterator1 rend () const {
01590                 return const_reverse_iterator1 (begin ());
01591             }
01592 #endif
01593 
01594             // Indices
01595             BOOST_UBLAS_INLINE
01596             size_type index1 () const {
01597                 return it2_.index1 ();
01598             }
01599             BOOST_UBLAS_INLINE
01600             size_type index2 () const {
01601                 return it2_.index2 ();
01602             }
01603 
01604             // Assignment
01605             BOOST_UBLAS_INLINE
01606             const_iterator2 &operator = (const const_iterator2 &it) {
01607                 container_const_reference<self_type>::assign (&it ());
01608                 it2_ = it.it2_;
01609                 return *this;
01610             }
01611 
01612             // Comparison
01613             BOOST_UBLAS_INLINE
01614             bool operator == (const const_iterator2 &it) const {
01615                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
01616                 return it2_ == it.it2_;
01617             }
01618             BOOST_UBLAS_INLINE
01619             bool operator < (const const_iterator2 &it) const {
01620                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
01621                 return it2_ < it.it2_;
01622             }
01623 
01624         private:
01625             const_subiterator2_type it2_;
01626         };
01627 #endif
01628 
01629         BOOST_UBLAS_INLINE
01630         const_iterator2 begin2 () const {
01631             return find2 (0, 0, 0);
01632         }
01633         BOOST_UBLAS_INLINE
01634         const_iterator2 end2 () const {
01635             return find2 (0, 0, size2 ());
01636         }
01637 
01638 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
01639         class iterator2:
01640             public container_reference<triangular_adaptor>,
01641             public random_access_iterator_base<typename iterator_restrict_traits<
01642                                                    typename subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,
01643                                                iterator2, value_type> {
01644         public:
01645             typedef typename subiterator2_type::value_type value_type;
01646             typedef typename subiterator2_type::difference_type difference_type;
01647             typedef typename subiterator2_type::reference reference;
01648             typedef typename subiterator2_type::pointer pointer;
01649 
01650             typedef iterator1 dual_iterator_type;
01651             typedef reverse_iterator1 dual_reverse_iterator_type;
01652 
01653             // Construction and destruction
01654             BOOST_UBLAS_INLINE
01655             iterator2 ():
01656                 container_reference<self_type> (), it2_ () {}
01657             BOOST_UBLAS_INLINE
01658             iterator2 (self_type &m, const subiterator2_type &it2):
01659                 container_reference<self_type> (m), it2_ (it2) {}
01660 
01661             // Arithmetic
01662             BOOST_UBLAS_INLINE
01663             iterator2 &operator ++ () {
01664                 ++ it2_;
01665                 return *this;
01666             }
01667             BOOST_UBLAS_INLINE
01668             iterator2 &operator -- () {
01669                 -- it2_;
01670                 return *this;
01671             }
01672             BOOST_UBLAS_INLINE
01673             iterator2 &operator += (difference_type n) {
01674                 it2_ += n;
01675                 return *this;
01676             }
01677             BOOST_UBLAS_INLINE
01678             iterator2 &operator -= (difference_type n) {
01679                 it2_ -= n;
01680                 return *this;
01681             }
01682             BOOST_UBLAS_INLINE
01683             difference_type operator - (const iterator2 &it) const {
01684                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
01685                 return it2_ - it.it2_;
01686             }
01687 
01688             // Dereference
01689             BOOST_UBLAS_INLINE
01690             reference operator * () const {
01691                 size_type i = index1 ();
01692                 size_type j = index2 ();
01693                 BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());
01694                 BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());
01695                 if (triangular_type::other (i, j))
01696                     return *it2_;
01697                 else
01698                     return (*this) () (i, j);
01699             }
01700             BOOST_UBLAS_INLINE
01701             reference operator [] (difference_type n) const {
01702                 return *(*this + n);
01703             }
01704 
01705 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
01706             BOOST_UBLAS_INLINE
01707 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
01708             typename self_type::
01709 #endif
01710             iterator1 begin () const {
01711                 return (*this) ().find1 (1, 0, index2 ());
01712             }
01713             BOOST_UBLAS_INLINE
01714 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
01715             typename self_type::
01716 #endif
01717             iterator1 end () const {
01718                 return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
01719             }
01720             BOOST_UBLAS_INLINE
01721 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
01722             typename self_type::
01723 #endif
01724             reverse_iterator1 rbegin () const {
01725                 return reverse_iterator1 (end ());
01726             }
01727             BOOST_UBLAS_INLINE
01728 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
01729             typename self_type::
01730 #endif
01731             reverse_iterator1 rend () const {
01732                 return reverse_iterator1 (begin ());
01733             }
01734 #endif
01735 
01736             // Indices
01737             BOOST_UBLAS_INLINE
01738             size_type index1 () const {
01739                 return it2_.index1 ();
01740             }
01741             BOOST_UBLAS_INLINE
01742             size_type index2 () const {
01743                 return it2_.index2 ();
01744             }
01745 
01746             // Assignment
01747             BOOST_UBLAS_INLINE
01748             iterator2 &operator = (const iterator2 &it) {
01749                 container_reference<self_type>::assign (&it ());
01750                 it2_ = it.it2_;
01751                 return *this;
01752             }
01753 
01754             // Comparison
01755             BOOST_UBLAS_INLINE
01756             bool operator == (const iterator2 &it) const {
01757                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
01758                 return it2_ == it.it2_;
01759             }
01760             BOOST_UBLAS_INLINE
01761             bool operator < (const iterator2 &it) const {
01762                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
01763                 return it2_ < it.it2_;
01764             }
01765 
01766         private:
01767             subiterator2_type it2_;
01768 
01769             friend class const_iterator2;
01770         };
01771 #endif
01772 
01773         BOOST_UBLAS_INLINE
01774         iterator2 begin2 () {
01775             return find2 (0, 0, 0);
01776         }
01777         BOOST_UBLAS_INLINE
01778         iterator2 end2 () {
01779             return find2 (0, 0, size2 ());
01780         }
01781 
01782         // Reverse iterators
01783 
01784         BOOST_UBLAS_INLINE
01785         const_reverse_iterator1 rbegin1 () const {
01786             return const_reverse_iterator1 (end1 ());
01787         }
01788         BOOST_UBLAS_INLINE
01789         const_reverse_iterator1 rend1 () const {
01790             return const_reverse_iterator1 (begin1 ());
01791         }
01792 
01793         BOOST_UBLAS_INLINE
01794         reverse_iterator1 rbegin1 () {
01795             return reverse_iterator1 (end1 ());
01796         }
01797         BOOST_UBLAS_INLINE
01798         reverse_iterator1 rend1 () {
01799             return reverse_iterator1 (begin1 ());
01800         }
01801 
01802         BOOST_UBLAS_INLINE
01803         const_reverse_iterator2 rbegin2 () const {
01804             return const_reverse_iterator2 (end2 ());
01805         }
01806         BOOST_UBLAS_INLINE
01807         const_reverse_iterator2 rend2 () const {
01808             return const_reverse_iterator2 (begin2 ());
01809         }
01810 
01811         BOOST_UBLAS_INLINE
01812         reverse_iterator2 rbegin2 () {
01813             return reverse_iterator2 (end2 ());
01814         }
01815         BOOST_UBLAS_INLINE
01816         reverse_iterator2 rend2 () {
01817             return reverse_iterator2 (begin2 ());
01818         }
01819 
01820     private:
01821         matrix_closure_type data_;
01822         static const value_type zero_;
01823         static const value_type one_;
01824     };
01825 
01826     template<class M, class TRI>
01827     const typename triangular_adaptor<M, TRI>::value_type triangular_adaptor<M, TRI>::zero_ = value_type/*zero*/();
01828     template<class M, class TRI>
01829     const typename triangular_adaptor<M, TRI>::value_type triangular_adaptor<M, TRI>::one_ (1);
01830 
01831     template <class M, class TRI>
01832     struct vector_temporary_traits< triangular_adaptor<M, TRI> >
01833     : vector_temporary_traits< typename boost::remove_const<M>::type > {} ;
01834     template <class M, class TRI>
01835     struct vector_temporary_traits< const triangular_adaptor<M, TRI> >
01836     : vector_temporary_traits< typename boost::remove_const<M>::type > {} ;
01837 
01838     template <class M, class TRI>
01839     struct matrix_temporary_traits< triangular_adaptor<M, TRI> >
01840     : matrix_temporary_traits< typename boost::remove_const<M>::type > {};
01841     template <class M, class TRI>
01842     struct matrix_temporary_traits< const triangular_adaptor<M, TRI> >
01843     : matrix_temporary_traits< typename boost::remove_const<M>::type > {};
01844 
01845 
01846     template<class E1, class E2>
01847     struct matrix_vector_solve_traits {
01848         typedef typename promote_traits<typename E1::value_type, typename E2::value_type>::promote_type promote_type;
01849         typedef vector<promote_type> result_type;
01850     };
01851 
01852     // Operations:
01853     //  n * (n - 1) / 2 + n = n * (n + 1) / 2 multiplications,
01854     //  n * (n - 1) / 2 additions
01855 
01856     // Dense (proxy) case
01857     template<class E1, class E2>
01858     BOOST_UBLAS_INLINE
01859     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
01860                         lower_tag, column_major_tag, dense_proxy_tag) {
01861         typedef typename E2::size_type size_type;
01862         typedef typename E2::difference_type difference_type;
01863         typedef typename E2::value_type value_type;
01864 
01865         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
01866         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
01867         size_type size = e2 ().size ();
01868         for (size_type n = 0; n < size; ++ n) {
01869 #ifndef BOOST_UBLAS_SINGULAR_CHECK
01870             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
01871 #else
01872             if (e1 () (n, n) == value_type/*zero*/())
01873                 singular ().raise ();
01874 #endif
01875             value_type t = e2 () (n) /= e1 () (n, n);
01876             if (t != value_type/*zero*/()) {
01877                 for (size_type m = n + 1; m < size; ++ m)
01878                     e2 () (m) -= e1 () (m, n) * t;
01879             }
01880         }
01881     }
01882     // Packed (proxy) case
01883     template<class E1, class E2>
01884     BOOST_UBLAS_INLINE
01885     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
01886                         lower_tag, column_major_tag, packed_proxy_tag) {
01887         typedef typename E2::size_type size_type;
01888         typedef typename E2::difference_type difference_type;
01889         typedef typename E2::value_type value_type;
01890 
01891         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
01892         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
01893         size_type size = e2 ().size ();
01894         for (size_type n = 0; n < size; ++ n) {
01895 #ifndef BOOST_UBLAS_SINGULAR_CHECK
01896             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
01897 #else
01898             if (e1 () (n, n) == value_type/*zero*/())
01899                 singular ().raise ();
01900 #endif
01901             value_type t = e2 () (n) /= e1 () (n, n);
01902             if (t != value_type/*zero*/()) {
01903                 typename E1::const_iterator1 it1e1 (e1 ().find1 (1, n + 1, n));
01904                 typename E1::const_iterator1 it1e1_end (e1 ().find1 (1, e1 ().size1 (), n));
01905                 difference_type m (it1e1_end - it1e1);
01906                 while (-- m >= 0)
01907                     e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1;
01908             }
01909         }
01910     }
01911     // Sparse (proxy) case
01912     template<class E1, class E2>
01913     BOOST_UBLAS_INLINE
01914     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
01915                         lower_tag, column_major_tag, unknown_storage_tag) {
01916         typedef typename E2::size_type size_type;
01917         typedef typename E2::difference_type difference_type;
01918         typedef typename E2::value_type value_type;
01919 
01920         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
01921         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
01922         size_type size = e2 ().size ();
01923         for (size_type n = 0; n < size; ++ n) {
01924 #ifndef BOOST_UBLAS_SINGULAR_CHECK
01925             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
01926 #else
01927             if (e1 () (n, n) == value_type/*zero*/())
01928                 singular ().raise ();
01929 #endif
01930             value_type t = e2 () (n) /= e1 () (n, n);
01931             if (t != value_type/*zero*/()) {
01932                 typename E1::const_iterator1 it1e1 (e1 ().find1 (1, n + 1, n));
01933                 typename E1::const_iterator1 it1e1_end (e1 ().find1 (1, e1 ().size1 (), n));
01934                 while (it1e1 != it1e1_end)
01935                     e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1;
01936             }
01937         }
01938     }
01939     // Redirectors :-)
01940     template<class E1, class E2>
01941     BOOST_UBLAS_INLINE
01942     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
01943                         lower_tag, column_major_tag) {
01944         typedef typename E1::storage_category storage_category;
01945         inplace_solve (e1, e2,
01946                        lower_tag (), column_major_tag (), storage_category ());
01947     }
01948     template<class E1, class E2>
01949     BOOST_UBLAS_INLINE
01950     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
01951                         lower_tag, row_major_tag) {
01952         typedef typename E1::storage_category storage_category;
01953         inplace_solve (e2, trans (e1),
01954                        upper_tag (), row_major_tag (), storage_category ());
01955     }
01956     // Dispatcher
01957     template<class E1, class E2>
01958     BOOST_UBLAS_INLINE
01959     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
01960                         lower_tag) {
01961         typedef typename E1::orientation_category orientation_category;
01962         inplace_solve (e1, e2,
01963                        lower_tag (), orientation_category ());
01964     }
01965     template<class E1, class E2>
01966     BOOST_UBLAS_INLINE
01967     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
01968                         unit_lower_tag) {
01969         typedef typename E1::orientation_category orientation_category;
01970         inplace_solve (triangular_adaptor<const E1, unit_lower> (e1 ()), e2,
01971                        unit_lower_tag (), orientation_category ());
01972     }
01973 
01974     // Dense (proxy) case
01975     template<class E1, class E2>
01976     BOOST_UBLAS_INLINE
01977     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
01978                         upper_tag, column_major_tag, dense_proxy_tag) {
01979         typedef typename E2::size_type size_type;
01980         typedef typename E2::difference_type difference_type;
01981         typedef typename E2::value_type value_type;
01982 
01983         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
01984         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
01985         size_type size = e2 ().size ();
01986         for (difference_type n = size - 1; n >= 0; -- n) {
01987 #ifndef BOOST_UBLAS_SINGULAR_CHECK
01988             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
01989 #else
01990             if (e1 () (n, n) == value_type/*zero*/())
01991                 singular ().raise ();
01992 #endif
01993             value_type t = e2 () (n) /= e1 () (n, n);
01994             if (t != value_type/*zero*/()) {
01995                 for (difference_type m = n - 1; m >= 0; -- m)
01996                     e2 () (m) -= e1 () (m, n) * t;
01997             }
01998         }
01999     }
02000     // Packed (proxy) case
02001     template<class E1, class E2>
02002     BOOST_UBLAS_INLINE
02003     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
02004                         upper_tag, column_major_tag, packed_proxy_tag) {
02005         typedef typename E2::size_type size_type;
02006         typedef typename E2::difference_type difference_type;
02007         typedef typename E2::value_type value_type;
02008 
02009         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
02010         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
02011         size_type size = e2 ().size ();
02012         for (difference_type n = size - 1; n >= 0; -- n) {
02013 #ifndef BOOST_UBLAS_SINGULAR_CHECK
02014             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
02015 #else
02016             if (e1 () (n, n) == value_type/*zero*/())
02017                 singular ().raise ();
02018 #endif
02019             value_type t = e2 () (n) /= e1 () (n, n);
02020             if (t != value_type/*zero*/()) {
02021                 typename E1::const_reverse_iterator1 it1e1 (e1 ().find1 (1, n, n));
02022                 typename E1::const_reverse_iterator1 it1e1_rend (e1 ().find1 (1, 0, n));
02023                 difference_type m (it1e1_rend - it1e1);
02024                 while (-- m >= 0)
02025                     e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1;
02026             }
02027         }
02028     }
02029     // Sparse (proxy) case
02030     template<class E1, class E2>
02031     BOOST_UBLAS_INLINE
02032     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
02033                         upper_tag, column_major_tag, unknown_storage_tag) {
02034         typedef typename E2::size_type size_type;
02035         typedef typename E2::difference_type difference_type;
02036         typedef typename E2::value_type value_type;
02037 
02038         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
02039         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
02040         size_type size = e2 ().size ();
02041         for (difference_type n = size - 1; n >= 0; -- n) {
02042 #ifndef BOOST_UBLAS_SINGULAR_CHECK
02043             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
02044 #else
02045             if (e1 () (n, n) == value_type/*zero*/())
02046                 singular ().raise ();
02047 #endif
02048             value_type t = e2 () (n) /= e1 () (n, n);
02049             if (t != value_type/*zero*/()) {
02050                 typename E1::const_reverse_iterator1 it1e1 (e1 ().find1 (1, n, n));
02051                 typename E1::const_reverse_iterator1 it1e1_rend (e1 ().find1 (1, 0, n));
02052                 while (it1e1 != it1e1_rend)
02053                     e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1;
02054             }
02055         }
02056     }
02057     // Redirectors :-)
02058     template<class E1, class E2>
02059     BOOST_UBLAS_INLINE
02060     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
02061                         upper_tag, column_major_tag) {
02062         typedef typename E1::storage_category storage_category;
02063         inplace_solve (e1, e2,
02064                        upper_tag (), column_major_tag (), storage_category ());
02065     }
02066     template<class E1, class E2>
02067     BOOST_UBLAS_INLINE
02068     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
02069                         upper_tag, row_major_tag) {
02070         typedef typename E1::storage_category storage_category;
02071         inplace_solve (e2, trans (e1),
02072                        lower_tag (), row_major_tag (), storage_category ());
02073     }
02074     // Dispatcher
02075     template<class E1, class E2>
02076     BOOST_UBLAS_INLINE
02077     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
02078                         upper_tag) {
02079         typedef typename E1::orientation_category orientation_category;
02080         inplace_solve (e1, e2,
02081                        upper_tag (), orientation_category ());
02082     }
02083     template<class E1, class E2>
02084     BOOST_UBLAS_INLINE
02085     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
02086                         unit_upper_tag) {
02087         typedef typename E1::orientation_category orientation_category;
02088         inplace_solve (triangular_adaptor<const E1, unit_upper> (e1 ()), e2,
02089                        unit_upper_tag (), orientation_category ());
02090     }
02091 
02092     template<class E1, class E2, class C>
02093     BOOST_UBLAS_INLINE
02094     typename matrix_vector_solve_traits<E1, E2>::result_type
02095     solve (const matrix_expression<E1> &e1,
02096            const vector_expression<E2> &e2,
02097            C) {
02098         typename matrix_vector_solve_traits<E1, E2>::result_type r (e2);
02099         inplace_solve (e1, r, C ());
02100         return r;
02101     }
02102 
02103     // Dense (proxy) case
02104     template<class E1, class E2>
02105     BOOST_UBLAS_INLINE
02106     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
02107                         lower_tag, row_major_tag, dense_proxy_tag) {
02108         typedef typename E1::size_type size_type;
02109         typedef typename E1::difference_type difference_type;
02110         typedef typename E1::value_type value_type;
02111 
02112         BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());
02113         BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());
02114         size_type size = e1 ().size ();
02115         for (difference_type n = size - 1; n >= 0; -- n) {
02116 #ifndef BOOST_UBLAS_SINGULAR_CHECK
02117             BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());
02118 #else
02119             if (e2 () (n, n) == value_type/*zero*/())
02120                 singular ().raise ();
02121 #endif
02122             value_type t = e1 () (n) /= e2 () (n, n);
02123             if (t != value_type/*zero*/()) {
02124                 for (difference_type m = n - 1; m >= 0; -- m)
02125                     e1 () (m) -= t * e2 () (n, m);
02126             }
02127         }
02128     }
02129     // Packed (proxy) case
02130     template<class E1, class E2>
02131     BOOST_UBLAS_INLINE
02132     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
02133                         lower_tag, row_major_tag, packed_proxy_tag) {
02134         typedef typename E1::size_type size_type;
02135         typedef typename E1::difference_type difference_type;
02136         typedef typename E1::value_type value_type;
02137 
02138         BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());
02139         BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());
02140         size_type size = e1 ().size ();
02141         for (difference_type n = size - 1; n >= 0; -- n) {
02142 #ifndef BOOST_UBLAS_SINGULAR_CHECK
02143             BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());
02144 #else
02145             if (e2 () (n, n) == value_type/*zero*/())
02146                 singular ().raise ();
02147 #endif
02148             value_type t = e1 () (n) /= e2 () (n, n);
02149             if (t != value_type/*zero*/()) {
02150                 typename E2::const_reverse_iterator2 it2e2 (e2 ().find2 (1, n, n));
02151                 typename E2::const_reverse_iterator2 it2e2_rend (e2 ().find2 (1, n, 0));
02152                 difference_type m (it2e2_rend - it2e2);
02153                 while (-- m >= 0)
02154                     e1 () (it2e2.index2 ()) -= *it2e2 * t, ++ it2e2;
02155             }
02156         }
02157     }
02158     // Sparse (proxy) case
02159     template<class E1, class E2>
02160     BOOST_UBLAS_INLINE
02161     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
02162                         lower_tag, row_major_tag, unknown_storage_tag) {
02163         typedef typename E1::size_type size_type;
02164         typedef typename E1::difference_type difference_type;
02165         typedef typename E1::value_type value_type;
02166 
02167         BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());
02168         BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());
02169         size_type size = e1 ().size ();
02170         for (difference_type n = size - 1; n >= 0; -- n) {
02171 #ifndef BOOST_UBLAS_SINGULAR_CHECK
02172             BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());
02173 #else
02174             if (e2 () (n, n) == value_type/*zero*/())
02175                 singular ().raise ();
02176 #endif
02177             value_type t = e1 () (n) /= e2 () (n, n);
02178             if (t != value_type/*zero*/()) {
02179                 typename E2::const_reverse_iterator2 it2e2 (e2 ().find2 (1, n, n));
02180                 typename E2::const_reverse_iterator2 it2e2_rend (e2 ().find2 (1, n, 0));
02181                 while (it2e2 != it2e2_rend)
02182                     e1 () (it2e2.index2 ()) -= *it2e2 * t, ++ it2e2;
02183             }
02184         }
02185     }
02186     // Redirectors :-)
02187     template<class E1, class E2>
02188     BOOST_UBLAS_INLINE
02189     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
02190                         lower_tag, row_major_tag) {
02191         typedef typename E1::storage_category storage_category;
02192         inplace_solve (e1, e2,
02193                        lower_tag (), row_major_tag (), storage_category ());
02194     }
02195     template<class E1, class E2>
02196     BOOST_UBLAS_INLINE
02197     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
02198                         lower_tag, column_major_tag) {
02199         typedef typename E1::storage_category storage_category;
02200         inplace_solve (trans (e2), e1,
02201                        upper_tag (), row_major_tag (), storage_category ());
02202     }
02203     // Dispatcher
02204     template<class E1, class E2>
02205     BOOST_UBLAS_INLINE
02206     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
02207                         lower_tag) {
02208         typedef typename E2::orientation_category orientation_category;
02209         inplace_solve (e1, e2,
02210                        lower_tag (), orientation_category ());
02211     }
02212     template<class E1, class E2>
02213     BOOST_UBLAS_INLINE
02214     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
02215                         unit_lower_tag) {
02216         typedef typename E2::orientation_category orientation_category;
02217         inplace_solve (e1, triangular_adaptor<const E2, unit_lower> (e2 ()),
02218                        unit_lower_tag (), orientation_category ());
02219     }
02220 
02221     // Dense (proxy) case
02222     template<class E1, class E2>
02223     BOOST_UBLAS_INLINE
02224     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
02225                         upper_tag, row_major_tag, dense_proxy_tag) {
02226         typedef typename E1::size_type size_type;
02227         typedef typename E1::difference_type difference_type;
02228         typedef typename E1::value_type value_type;
02229 
02230         BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());
02231         BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());
02232         size_type size = e1 ().size ();
02233         for (size_type n = 0; n < size; ++ n) {
02234 #ifndef BOOST_UBLAS_SINGULAR_CHECK
02235             BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());
02236 #else
02237             if (e2 () (n, n) == value_type/*zero*/())
02238                 singular ().raise ();
02239 #endif
02240             value_type t = e1 () (n) /= e2 () (n, n);
02241             if (t != value_type/*zero*/()) {
02242                 for (size_type m = n + 1; m < size; ++ m)
02243                     e1 () (m) -= t * e2 () (n, m);
02244             }
02245         }
02246     }
02247     // Packed (proxy) case
02248     template<class E1, class E2>
02249     BOOST_UBLAS_INLINE
02250     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
02251                         upper_tag, row_major_tag, packed_proxy_tag) {
02252         typedef typename E1::size_type size_type;
02253         typedef typename E1::difference_type difference_type;
02254         typedef typename E1::value_type value_type;
02255 
02256         BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());
02257         BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());
02258         size_type size = e1 ().size ();
02259         for (size_type n = 0; n < size; ++ n) {
02260 #ifndef BOOST_UBLAS_SINGULAR_CHECK
02261             BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());
02262 #else
02263             if (e2 () (n, n) == value_type/*zero*/())
02264                 singular ().raise ();
02265 #endif
02266             value_type t = e1 () (n) /= e2 () (n, n);
02267             if (t != value_type/*zero*/()) {
02268                 typename E2::const_iterator2 it2e2 (e2 ().find2 (1, n, n + 1));
02269                 typename E2::const_iterator2 it2e2_end (e2 ().find2 (1, n, e2 ().size2 ()));
02270                 difference_type m (it2e2_end - it2e2);
02271                 while (-- m >= 0)
02272                     e1 () (it2e2.index2 ()) -= *it2e2 * t, ++ it2e2;
02273             }
02274         }
02275     }
02276     // Sparse (proxy) case
02277     template<class E1, class E2>
02278     BOOST_UBLAS_INLINE
02279     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
02280                         upper_tag, row_major_tag, unknown_storage_tag) {
02281         typedef typename E1::size_type size_type;
02282         typedef typename E1::difference_type difference_type;
02283         typedef typename E1::value_type value_type;
02284 
02285         BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());
02286         BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());
02287         size_type size = e1 ().size ();
02288         for (size_type n = 0; n < size; ++ n) {
02289 #ifndef BOOST_UBLAS_SINGULAR_CHECK
02290             BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());
02291 #else
02292             if (e2 () (n, n) == value_type/*zero*/())
02293                 singular ().raise ();
02294 #endif
02295             value_type t = e1 () (n) /= e2 () (n, n);
02296             if (t != value_type/*zero*/()) {
02297                 typename E2::const_iterator2 it2e2 (e2 ().find2 (1, n, n + 1));
02298                 typename E2::const_iterator2 it2e2_end (e2 ().find2 (1, n, e2 ().size2 ()));
02299                 while (it2e2 != it2e2_end)
02300                     e1 () (it2e2.index2 ()) -= *it2e2 * t, ++ it2e2;
02301             }
02302         }
02303     }
02304     // Redirectors :-)
02305     template<class E1, class E2>
02306     BOOST_UBLAS_INLINE
02307     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
02308                         upper_tag, row_major_tag) {
02309         typedef typename E1::storage_category storage_category;
02310         inplace_solve (e1, e2,
02311                        upper_tag (), row_major_tag (), storage_category ());
02312     }
02313     template<class E1, class E2>
02314     BOOST_UBLAS_INLINE
02315     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
02316                         upper_tag, column_major_tag) {
02317         typedef typename E1::storage_category storage_category;
02318         inplace_solve (trans (e2), e1,
02319                        lower_tag (), row_major_tag (), storage_category ());
02320     }
02321     // Dispatcher
02322     template<class E1, class E2>
02323     BOOST_UBLAS_INLINE
02324     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
02325                         upper_tag) {
02326         typedef typename E2::orientation_category orientation_category;
02327         inplace_solve (e1, e2,
02328                        upper_tag (), orientation_category ());
02329     }
02330     template<class E1, class E2>
02331     BOOST_UBLAS_INLINE
02332     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
02333                         unit_upper_tag) {
02334         typedef typename E2::orientation_category orientation_category;
02335         inplace_solve (e1, triangular_adaptor<const E2, unit_upper> (e2 ()),
02336                        unit_upper_tag (), orientation_category ());
02337     }
02338 
02339     template<class E1, class E2, class C>
02340     BOOST_UBLAS_INLINE
02341     typename matrix_vector_solve_traits<E1, E2>::result_type
02342     solve (const vector_expression<E1> &e1,
02343            const matrix_expression<E2> &e2,
02344            C) {
02345         typename matrix_vector_solve_traits<E1, E2>::result_type r (e1);
02346         inplace_solve (r, e2, C ());
02347         return r;
02348     }
02349 
02350     template<class E1, class E2>
02351     struct matrix_matrix_solve_traits {
02352         typedef typename promote_traits<typename E1::value_type, typename E2::value_type>::promote_type promote_type;
02353         typedef matrix<promote_type> result_type;
02354     };
02355 
02356     // Operations:
02357     //  k * n * (n - 1) / 2 + k * n = k * n * (n + 1) / 2 multiplications,
02358     //  k * n * (n - 1) / 2 additions
02359 
02360     // Dense (proxy) case
02361     template<class E1, class E2>
02362     BOOST_UBLAS_INLINE
02363     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
02364                         lower_tag, dense_proxy_tag) {
02365         typedef typename E2::size_type size_type;
02366         typedef typename E2::difference_type difference_type;
02367         typedef typename E2::value_type value_type;
02368 
02369         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
02370         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());
02371         size_type size1 = e2 ().size1 ();
02372         size_type size2 = e2 ().size2 ();
02373         for (size_type n = 0; n < size1; ++ n) {
02374 #ifndef BOOST_UBLAS_SINGULAR_CHECK
02375             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
02376 #else
02377             if (e1 () (n, n) == value_type/*zero*/())
02378                 singular ().raise ();
02379 #endif
02380             for (size_type l = 0; l < size2; ++ l) {
02381                 value_type t = e2 () (n, l) /= e1 () (n, n);
02382                 if (t != value_type/*zero*/()) {
02383                     for (size_type m = n + 1; m < size1; ++ m)
02384                         e2 () (m, l) -= e1 () (m, n) * t;
02385                 }
02386             }
02387         }
02388     }
02389     // Packed (proxy) case
02390     template<class E1, class E2>
02391     BOOST_UBLAS_INLINE
02392     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
02393                         lower_tag, packed_proxy_tag) {
02394         typedef typename E2::size_type size_type;
02395         typedef typename E2::difference_type difference_type;
02396         typedef typename E2::value_type value_type;
02397 
02398         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
02399         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());
02400         size_type size1 = e2 ().size1 ();
02401         size_type size2 = e2 ().size2 ();
02402         for (size_type n = 0; n < size1; ++ n) {
02403 #ifndef BOOST_UBLAS_SINGULAR_CHECK
02404             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
02405 #else
02406             if (e1 () (n, n) == value_type/*zero*/())
02407                 singular ().raise ();
02408 #endif
02409             for (size_type l = 0; l < size2; ++ l) {
02410                 value_type t = e2 () (n, l) /= e1 () (n, n);
02411                 if (t != value_type/*zero*/()) {
02412                     typename E1::const_iterator1 it1e1 (e1 ().find1 (1, n + 1, n));
02413                     typename E1::const_iterator1 it1e1_end (e1 ().find1 (1, e1 ().size1 (), n));
02414                     difference_type m (it1e1_end - it1e1);
02415                     while (-- m >= 0)
02416                         e2 () (it1e1.index1 (), l) -= *it1e1 * t, ++ it1e1;
02417                 }
02418             }
02419         }
02420     }
02421     // Sparse (proxy) case
02422     template<class E1, class E2>
02423     BOOST_UBLAS_INLINE
02424     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
02425                         lower_tag, unknown_storage_tag) {
02426         typedef typename E2::size_type size_type;
02427         typedef typename E2::difference_type difference_type;
02428         typedef typename E2::value_type value_type;
02429 
02430         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
02431         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());
02432         size_type size1 = e2 ().size1 ();
02433         size_type size2 = e2 ().size2 ();
02434         for (size_type n = 0; n < size1; ++ n) {
02435 #ifndef BOOST_UBLAS_SINGULAR_CHECK
02436             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
02437 #else
02438             if (e1 () (n, n) == value_type/*zero*/())
02439                 singular ().raise ();
02440 #endif
02441             for (size_type l = 0; l < size2; ++ l) {
02442                 value_type t = e2 () (n, l) /= e1 () (n, n);
02443                 if (t != value_type/*zero*/()) {
02444                     typename E1::const_iterator1 it1e1 (e1 ().find1 (1, n + 1, n));
02445                     typename E1::const_iterator1 it1e1_end (e1 ().find1 (1, e1 ().size1 (), n));
02446                     while (it1e1 != it1e1_end)
02447                         e2 () (it1e1.index1 (), l) -= *it1e1 * t, ++ it1e1;
02448                 }
02449             }
02450         }
02451     }
02452     // Dispatcher
02453     template<class E1, class E2>
02454     BOOST_UBLAS_INLINE
02455     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
02456                         lower_tag) {
02457         typedef typename E1::storage_category dispatch_category;
02458         inplace_solve (e1, e2,
02459                        lower_tag (), dispatch_category ());
02460     }
02461     template<class E1, class E2>
02462     BOOST_UBLAS_INLINE
02463     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
02464                         unit_lower_tag) {
02465         typedef typename E1::storage_category dispatch_category;
02466         inplace_solve (triangular_adaptor<const E1, unit_lower> (e1 ()), e2,
02467                        unit_lower_tag (), dispatch_category ());
02468     }
02469 
02470     // Dense (proxy) case
02471     template<class E1, class E2>
02472     BOOST_UBLAS_INLINE
02473     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
02474                         upper_tag, dense_proxy_tag) {
02475         typedef typename E2::size_type size_type;
02476         typedef typename E2::difference_type difference_type;
02477         typedef typename E2::value_type value_type;
02478 
02479         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
02480         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());
02481         size_type size1 = e2 ().size1 ();
02482         size_type size2 = e2 ().size2 ();
02483         for (difference_type n = size1 - 1; n >= 0; -- n) {
02484 #ifndef BOOST_UBLAS_SINGULAR_CHECK
02485             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
02486 #else
02487             if (e1 () (n, n) == value_type/*zero*/())
02488                 singular ().raise ();
02489 #endif
02490             for (difference_type l = size2 - 1; l >= 0; -- l) {
02491                 value_type t = e2 () (n, l) /= e1 () (n, n);
02492                 if (t != value_type/*zero*/()) {
02493                     for (difference_type m = n - 1; m >= 0; -- m)
02494                         e2 () (m, l) -= e1 () (m, n) * t;
02495                 }
02496             }
02497         }
02498     }
02499     // Packed (proxy) case
02500     template<class E1, class E2>
02501     BOOST_UBLAS_INLINE
02502     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
02503                         upper_tag, packed_proxy_tag) {
02504         typedef typename E2::size_type size_type;
02505         typedef typename E2::difference_type difference_type;
02506         typedef typename E2::value_type value_type;
02507 
02508         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
02509         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());
02510         size_type size1 = e2 ().size1 ();
02511         size_type size2 = e2 ().size2 ();
02512         for (difference_type n = size1 - 1; n >= 0; -- n) {
02513 #ifndef BOOST_UBLAS_SINGULAR_CHECK
02514             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
02515 #else
02516             if (e1 () (n, n) == value_type/*zero*/())
02517                 singular ().raise ();
02518 #endif
02519             for (difference_type l = size2 - 1; l >= 0; -- l) {
02520                 value_type t = e2 () (n, l) /= e1 () (n, n);
02521                 if (t != value_type/*zero*/()) {
02522                     typename E1::const_reverse_iterator1 it1e1 (e1 ().find1 (1, n, n));
02523                     typename E1::const_reverse_iterator1 it1e1_rend (e1 ().find1 (1, 0, n));
02524                     difference_type m (it1e1_rend - it1e1);
02525                     while (-- m >= 0)
02526                         e2 () (it1e1.index1 (), l) -= *it1e1 * t, ++ it1e1;
02527                 }
02528             }
02529         }
02530     }
02531     // Sparse (proxy) case
02532     template<class E1, class E2>
02533     BOOST_UBLAS_INLINE
02534     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
02535                         upper_tag, unknown_storage_tag) {
02536         typedef typename E2::size_type size_type;
02537         typedef typename E2::difference_type difference_type;
02538         typedef typename E2::value_type value_type;
02539 
02540         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
02541         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());
02542         size_type size1 = e2 ().size1 ();
02543         size_type size2 = e2 ().size2 ();
02544         for (difference_type n = size1 - 1; n >= 0; -- n) {
02545 #ifndef BOOST_UBLAS_SINGULAR_CHECK
02546             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
02547 #else
02548             if (e1 () (n, n) == value_type/*zero*/())
02549                 singular ().raise ();
02550 #endif
02551             for (difference_type l = size2 - 1; l >= 0; -- l) {
02552                 value_type t = e2 () (n, l) /= e1 () (n, n);
02553                 if (t != value_type/*zero*/()) {
02554                     typename E1::const_reverse_iterator1 it1e1 (e1 ().find1 (1, n, n));
02555                     typename E1::const_reverse_iterator1 it1e1_rend (e1 ().find1 (1, 0, n));
02556                     while (it1e1 != it1e1_rend)
02557                         e2 () (it1e1.index1 (), l) -= *it1e1 * t, ++ it1e1;
02558                 }
02559             }
02560         }
02561     }
02562     // Dispatcher
02563     template<class E1, class E2>
02564     BOOST_UBLAS_INLINE
02565     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
02566                         upper_tag) {
02567         typedef typename E1::storage_category dispatch_category;
02568         inplace_solve (e1, e2,
02569                        upper_tag (), dispatch_category ());
02570     }
02571     template<class E1, class E2>
02572     BOOST_UBLAS_INLINE
02573     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
02574                         unit_upper_tag) {
02575         typedef typename E1::storage_category dispatch_category;
02576         inplace_solve (triangular_adaptor<const E1, unit_upper> (e1 ()), e2,
02577                        unit_upper_tag (), dispatch_category ());
02578     }
02579 
02580     template<class E1, class E2, class C>
02581     BOOST_UBLAS_INLINE
02582     typename matrix_matrix_solve_traits<E1, E2>::result_type
02583     solve (const matrix_expression<E1> &e1,
02584            const matrix_expression<E2> &e2,
02585            C) {
02586         typename matrix_matrix_solve_traits<E1, E2>::result_type r (e2);
02587         inplace_solve (e1, r, C ());
02588         return r;
02589     }
02590 
02591 }}}
02592 
02593 #endif