- Main Page
- Related Pages
- Classes
- Files

- File List

# triangular.hpp

    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