boost/numeric/ublas/triangular.hpp

   1 //
   2 //  Copyright (c) 2000-2002
   3 //  Joerg Walter, Mathias Koch
   4 //
   5 //  Distributed under the Boost Software License, Version 1.0. (See
   6 //  accompanying file LICENSE_1_0.txt or copy at
   7 //  http://www.boost.org/LICENSE_1_0.txt)
   8 //
   9 //  The authors gratefully acknowledge the support of
  10 //  GeNeSys mbH & Co. KG in producing this work.
  11 //
  12
  13 #ifndef _BOOST_UBLAS_TRIANGULAR_
  14 #define _BOOST_UBLAS_TRIANGULAR_
  15
  16 #include <boost/numeric/ublas/matrix.hpp>
  17 #include <boost/numeric/ublas/detail/temporary.hpp>
  18 #include <boost/type_traits/remove_const.hpp>
  19
  20 // Iterators based on ideas of Jeremy Siek
  21
  22 namespace boost { namespace numeric { namespace ublas {
  23
  24     namespace detail {
  25         using namespace boost::numeric::ublas;
  26
  27         // Matrix resizing algorithm
  28         template <class L, class T, class M>
  29         BOOST_UBLAS_INLINE
  30         void matrix_resize_preserve (M& m, M& temporary) {
  31             typedef L layout_type;
  32             typedef T triangular_type;
  33             typedef typename M::size_type size_type;
  34             const size_type msize1 (m.size1 ());        // original size
  35             const size_type msize2 (m.size2 ());
  36             const size_type size1 (temporary.size1 ());    // new size is specified by temporary
  37             const size_type size2 (temporary.size2 ());
  38             // Common elements to preserve
  39             const size_type size1_min = (std::min) (size1, msize1);
  40             const size_type size2_min = (std::min) (size2, msize2);
  41             // Order for major and minor sizes
  42             const size_type major_size = layout_type::size_M (size1_min, size2_min);
  43             const size_type minor_size = layout_type::size_m (size1_min, size2_min);
  44             // Indexing copy over major
  45             for (size_type major = 0; major != major_size; ++major) {
  46                 for (size_type minor = 0; minor != minor_size; ++minor) {
  47                         // find indexes - use invertability of element_ functions
  48                     const size_type i1 = layout_type::index_M(major, minor);
  49                     const size_type i2 = layout_type::index_m(major, minor);
  50                     if ( triangular_type::other(i1,i2) ) {
  51                         temporary.data () [triangular_type::element (layout_type (), i1, size1, i2, size2)] =
  52                             m.data() [triangular_type::element (layout_type (), i1, msize1, i2, msize2)];
  53                     }
  54                 }
  55             }
  56             m.assign_temporary (temporary);
  57         }
  58     }
  59
  60     /** \brief A triangular matrix of values of type \c T.
  61      *
  62      * For a \f$(n \times n )\f$-dimensional lower triangular matrix and if \f$0 \leq i < n\f$, \f$0 \leq j < n\f$ and \f$i>j\f$ holds,
  63      * \f$m_{i,j}=0\f$. Furthermore if \f$m_{i,i}=1\f$, the matrix is called unit lower triangular.
  64      *
  65      * For a \f$(n \times n )\f$-dimensional upper triangular matrix and if \f$0 \leq i < n\f$, \f$0 \leq j < n\f$ and \f$i<j\f$ holds,
  66      * \f$m_{i,j}=0\f$. Furthermore if \f$m_{i,i}=1\f$, the matrix is called unit upper triangular.
  67      *
  68      * The default storage for triangular matrices is packed. Orientation and storage can also be specified.
  69      * Default is \c row_major and and unbounded_array. It is \b not required by the storage to initialize
  70      * elements of the matrix.
  71      *
  72      * \tparam T the type of object stored in the matrix (like double, float, complex, etc...)
  73      * \tparam TRI the type of the triangular matrix. It can either be \c lower or \c upper. Default is \c lower
  74      * \tparam L the storage organization. It can be either \c row_major or \c column_major. Default is \c row_major
  75      * \tparam A the type of Storage array. Default is \c unbounded_array
  76      */
  77     template<class T, class TRI, class L, class A>
  78     class triangular_matrix:
  79         public matrix_container<triangular_matrix<T, TRI, L, A> > {
  80
  81         typedef T *pointer;
  82         typedef TRI triangular_type;
  83         typedef L layout_type;
  84         typedef triangular_matrix<T, TRI, L, A> self_type;
  85     public:
  86 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
  87         using matrix_container<self_type>::operator ();
  88 #endif
  89         typedef typename A::size_type size_type;
  90         typedef typename A::difference_type difference_type;
  91         typedef T value_type;
  92         typedef const T &const_reference;
  93         typedef T &reference;
  94         typedef A array_type;
  95
  96         typedef const matrix_reference<const self_type> const_closure_type;
  97         typedef matrix_reference<self_type> closure_type;
  98         typedef vector<T, A> vector_temporary_type;
  99         typedef matrix<T, L, A> matrix_temporary_type;  // general sub-matrix
 100         typedef packed_tag storage_category;
 101         typedef typename L::orientation_category orientation_category;
 102
 103         // Construction and destruction
 104         BOOST_UBLAS_INLINE
 105         triangular_matrix ():
 106             matrix_container<self_type> (),
 107             size1_ (0), size2_ (0), data_ (0) {}
 108         BOOST_UBLAS_INLINE
 109         triangular_matrix (size_type size1, size_type size2):
 110             matrix_container<self_type> (),
 111             size1_ (size1), size2_ (size2), data_ (triangular_type::packed_size (layout_type (), size1, size2)) {
 112         }
 113         BOOST_UBLAS_INLINE
 114         triangular_matrix (size_type size1, size_type size2, const array_type &data):
 115             matrix_container<self_type> (),
 116             size1_ (size1), size2_ (size2), data_ (data) {}
 117         BOOST_UBLAS_INLINE
 118         triangular_matrix (const triangular_matrix &m):
 119             matrix_container<self_type> (),
 120             size1_ (m.size1_), size2_ (m.size2_), data_ (m.data_) {}
 121         template<class AE>
 122         BOOST_UBLAS_INLINE
 123         triangular_matrix (const matrix_expression<AE> &ae):
 124             matrix_container<self_type> (),
 125             size1_ (ae ().size1 ()), size2_ (ae ().size2 ()),
 126             data_ (triangular_type::packed_size (layout_type (), size1_, size2_)) {
 127             matrix_assign<scalar_assign> (*this, ae);
 128         }
 129
 130         // Accessors
 131         BOOST_UBLAS_INLINE
 132         size_type size1 () const {
 133             return size1_;
 134         }
 135         BOOST_UBLAS_INLINE
 136         size_type size2 () const {
 137             return size2_;
 138         }
 139
 140         // Storage accessors
 141         BOOST_UBLAS_INLINE
 142         const array_type &data () const {
 143             return data_;
 144         }
 145         BOOST_UBLAS_INLINE
 146         array_type &data () {
 147             return data_;
 148         }
 149
 150         // Resizing
 151         BOOST_UBLAS_INLINE
 152         void resize (size_type size1, size_type size2, bool preserve = true) {
 153             if (preserve) {
 154                 self_type temporary (size1, size2);
 155                 detail::matrix_resize_preserve<layout_type, triangular_type> (*this, temporary);
 156             }
 157             else {
 158                 data ().resize (triangular_type::packed_size (layout_type (), size1, size2));
 159                 size1_ = size1;
 160                 size2_ = size2;
 161             }
 162         }
 163         BOOST_UBLAS_INLINE
 164         void resize_packed_preserve (size_type size1, size_type size2) {
 165             size1_ = size1;
 166             size2_ = size2;
 167             data ().resize (triangular_type::packed_size (layout_type (), size1_, size2_), value_type ());
 168         }
 169
 170         // Element access
 171         BOOST_UBLAS_INLINE
 172         const_reference operator () (size_type i, size_type j) const {
 173             BOOST_UBLAS_CHECK (i < size1_, bad_index ());
 174             BOOST_UBLAS_CHECK (j < size2_, bad_index ());
 175             if (triangular_type::other (i, j))
 176                 return data () [triangular_type::element (layout_type (), i, size1_, j, size2_)];
 177             else if (triangular_type::one (i, j))
 178                 return one_;
 179             else
 180                 return zero_;
 181         }
 182         BOOST_UBLAS_INLINE
 183         reference at_element (size_type i, size_type j) {
 184             BOOST_UBLAS_CHECK (i < size1_, bad_index ());
 185             BOOST_UBLAS_CHECK (j < size2_, bad_index ());
 186             return data () [triangular_type::element (layout_type (), i, size1_, j, size2_)];
 187         }
 188         BOOST_UBLAS_INLINE
 189         reference operator () (size_type i, size_type j) {
 190             BOOST_UBLAS_CHECK (i < size1_, bad_index ());
 191             BOOST_UBLAS_CHECK (j < size2_, bad_index ());
 192             if (!triangular_type::other (i, j)) {
 193                 bad_index ().raise ();
 194                 // NEVER reached
 195             }
 196             return data () [triangular_type::element (layout_type (), i, size1_, j, size2_)];
 197         }
 198
 199         // Element assignment
 200         BOOST_UBLAS_INLINE
 201         reference insert_element (size_type i, size_type j, const_reference t) {
 202             return (operator () (i, j) = t);
 203         }
 204         BOOST_UBLAS_INLINE
 205         void erase_element (size_type i, size_type j) {
 206             operator () (i, j) = value_type/*zero*/();
 207         }
 208
 209         // Zeroing
 210         BOOST_UBLAS_INLINE
 211         void clear () {
 212             // data ().clear ();
 213             std::fill (data ().begin (), data ().end (), value_type/*zero*/());
 214         }
 215
 216         // Assignment
 217         BOOST_UBLAS_INLINE
 218         triangular_matrix &operator = (const triangular_matrix &m) {
 219             size1_ = m.size1_;
 220             size2_ = m.size2_;
 221             data () = m.data ();
 222             return *this;
 223         }
 224         BOOST_UBLAS_INLINE
 225         triangular_matrix &assign_temporary (triangular_matrix &m) {
 226             swap (m);
 227             return *this;
 228         }
 229         template<class AE>
 230         BOOST_UBLAS_INLINE
 231         triangular_matrix &operator = (const matrix_expression<AE> &ae) {
 232             self_type temporary (ae);
 233             return assign_temporary (temporary);
 234         }
 235         template<class AE>
 236         BOOST_UBLAS_INLINE
 237         triangular_matrix &assign (const matrix_expression<AE> &ae) {
 238             matrix_assign<scalar_assign> (*this, ae);
 239             return *this;
 240         }
 241         template<class AE>
 242         BOOST_UBLAS_INLINE
 243         triangular_matrix& operator += (const matrix_expression<AE> &ae) {
 244             self_type temporary (*this + ae);
 245             return assign_temporary (temporary);
 246         }
 247         template<class AE>
 248         BOOST_UBLAS_INLINE
 249         triangular_matrix &plus_assign (const matrix_expression<AE> &ae) {
 250             matrix_assign<scalar_plus_assign> (*this, ae);
 251             return *this;
 252         }
 253         template<class AE>
 254         BOOST_UBLAS_INLINE
 255         triangular_matrix& operator -= (const matrix_expression<AE> &ae) {
 256             self_type temporary (*this - ae);
 257             return assign_temporary (temporary);
 258         }
 259         template<class AE>
 260         BOOST_UBLAS_INLINE
 261         triangular_matrix &minus_assign (const matrix_expression<AE> &ae) {
 262             matrix_assign<scalar_minus_assign> (*this, ae);
 263             return *this;
 264         }
 265         template<class AT>
 266         BOOST_UBLAS_INLINE
 267         triangular_matrix& operator *= (const AT &at) {
 268             matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
 269             return *this;
 270         }
 271         template<class AT>
 272         BOOST_UBLAS_INLINE
 273         triangular_matrix& operator /= (const AT &at) {
 274             matrix_assign_scalar<scalar_divides_assign> (*this, at);
 275             return *this;
 276         }
 277
 278         // Swapping
 279         BOOST_UBLAS_INLINE
 280         void swap (triangular_matrix &m) {
 281             if (this != &m) {
 282                 // BOOST_UBLAS_CHECK (size2_ == m.size2_, bad_size ());
 283                 std::swap (size1_, m.size1_);
 284                 std::swap (size2_, m.size2_);
 285                 data ().swap (m.data ());
 286             }
 287         }
 288         BOOST_UBLAS_INLINE
 289         friend void swap (triangular_matrix &m1, triangular_matrix &m2) {
 290             m1.swap (m2);
 291         }
 292
 293         // Iterator types
 294 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
 295         typedef indexed_iterator1<self_type, packed_random_access_iterator_tag> iterator1;
 296         typedef indexed_iterator2<self_type, packed_random_access_iterator_tag> iterator2;
 297         typedef indexed_const_iterator1<self_type, packed_random_access_iterator_tag> const_iterator1;
 298         typedef indexed_const_iterator2<self_type, packed_random_access_iterator_tag> const_iterator2;
 299 #else
 300         class const_iterator1;
 301         class iterator1;
 302         class const_iterator2;
 303         class iterator2;
 304 #endif
 305         typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
 306         typedef reverse_iterator_base1<iterator1> reverse_iterator1;
 307         typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
 308         typedef reverse_iterator_base2<iterator2> reverse_iterator2;
 309
 310         // Element lookup
 311         BOOST_UBLAS_INLINE
 312         const_iterator1 find1 (int rank, size_type i, size_type j) const {
 313             if (rank == 1)
 314                 i = triangular_type::restrict1 (i, j, size1_, size2_);
 315             if (rank == 0)
 316                 i = triangular_type::global_restrict1 (i, size1_, j, size2_);
 317             return const_iterator1 (*this, i, j);
 318         }
 319         BOOST_UBLAS_INLINE
 320         iterator1 find1 (int rank, size_type i, size_type j) {
 321             if (rank == 1)
 322                 i = triangular_type::mutable_restrict1 (i, j, size1_, size2_);
 323             if (rank == 0)
 324                 i = triangular_type::global_mutable_restrict1 (i, size1_, j, size2_);
 325             return iterator1 (*this, i, j);
 326         }
 327         BOOST_UBLAS_INLINE
 328         const_iterator2 find2 (int rank, size_type i, size_type j) const {
 329             if (rank == 1)
 330                 j = triangular_type::restrict2 (i, j, size1_, size2_);
 331             if (rank == 0)
 332                 j = triangular_type::global_restrict2 (i, size1_, j, size2_);
 333             return const_iterator2 (*this, i, j);
 334         }
 335         BOOST_UBLAS_INLINE
 336         iterator2 find2 (int rank, size_type i, size_type j) {
 337             if (rank == 1)
 338                 j = triangular_type::mutable_restrict2 (i, j, size1_, size2_);
 339             if (rank == 0)
 340                 j = triangular_type::global_mutable_restrict2 (i, size1_, j, size2_);
 341             return iterator2 (*this, i, j);
 342         }
 343
 344         // Iterators simply are indices.
 345
 346 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
 347         class const_iterator1:
 348             public container_const_reference<triangular_matrix>,
 349             public random_access_iterator_base<packed_random_access_iterator_tag,
 350                                                const_iterator1, value_type> {
 351         public:
 352             typedef typename triangular_matrix::value_type value_type;
 353             typedef typename triangular_matrix::difference_type difference_type;
 354             typedef typename triangular_matrix::const_reference reference;
 355             typedef const typename triangular_matrix::pointer pointer;
 356
 357             typedef const_iterator2 dual_iterator_type;
 358             typedef const_reverse_iterator2 dual_reverse_iterator_type;
 359
 360             // Construction and destruction
 361             BOOST_UBLAS_INLINE
 362             const_iterator1 ():
 363                 container_const_reference<self_type> (), it1_ (), it2_ () {}
 364             BOOST_UBLAS_INLINE
 365             const_iterator1 (const self_type &m, size_type it1, size_type it2):
 366                 container_const_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
 367             BOOST_UBLAS_INLINE
 368             const_iterator1 (const iterator1 &it):
 369                 container_const_reference<self_type> (it ()), it1_ (it.it1_), it2_ (it.it2_) {}
 370
 371             // Arithmetic
 372             BOOST_UBLAS_INLINE
 373             const_iterator1 &operator ++ () {
 374                 ++ it1_;
 375                 return *this;
 376             }
 377             BOOST_UBLAS_INLINE
 378             const_iterator1 &operator -- () {
 379                 -- it1_;
 380                 return *this;
 381             }
 382             BOOST_UBLAS_INLINE
 383             const_iterator1 &operator += (difference_type n) {
 384                 it1_ += n;
 385                 return *this;
 386             }
 387             BOOST_UBLAS_INLINE
 388             const_iterator1 &operator -= (difference_type n) {
 389                 it1_ -= n;
 390                 return *this;
 391             }
 392             BOOST_UBLAS_INLINE
 393             difference_type operator - (const const_iterator1 &it) const {
 394                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
 395                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
 396                 return it1_ - it.it1_;
 397             }
 398
 399             // Dereference
 400             BOOST_UBLAS_INLINE
 401             const_reference operator * () const {
 402                 return (*this) () (it1_, it2_);
 403             }
 404             BOOST_UBLAS_INLINE
 405             const_reference operator [] (difference_type n) const {
 406                 return *(*this + n);
 407             }
 408
 409 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
 410             BOOST_UBLAS_INLINE
 411 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 412             typename self_type::
 413 #endif
 414             const_iterator2 begin () const {
 415                 return (*this) ().find2 (1, it1_, 0);
 416             }
 417             BOOST_UBLAS_INLINE
 418 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 419             typename self_type::
 420 #endif
 421             const_iterator2 cbegin () const {
 422                 return begin ();
 423             }
 424             BOOST_UBLAS_INLINE
 425 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 426             typename self_type::
 427 #endif
 428             const_iterator2 end () const {
 429                 return (*this) ().find2 (1, it1_, (*this) ().size2 ());
 430             }
 431             BOOST_UBLAS_INLINE
 432 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 433             typename self_type::
 434 #endif
 435             const_iterator2 cend () const {
 436                 return end ();
 437             }
 438             BOOST_UBLAS_INLINE
 439 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 440             typename self_type::
 441 #endif
 442             const_reverse_iterator2 rbegin () const {
 443                 return const_reverse_iterator2 (end ());
 444             }
 445             BOOST_UBLAS_INLINE
 446 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 447             typename self_type::
 448 #endif
 449             const_reverse_iterator2 crbegin () const {
 450                 return rbegin ();
 451             }
 452             BOOST_UBLAS_INLINE
 453 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 454             typename self_type::
 455 #endif
 456             const_reverse_iterator2 rend () const {
 457                 return const_reverse_iterator2 (begin ());
 458             }
 459             BOOST_UBLAS_INLINE
 460 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 461             typename self_type::
 462 #endif
 463             const_reverse_iterator2 crend () const {
 464                 return rend ();
 465             }
 466 #endif
 467
 468             // Indices
 469             BOOST_UBLAS_INLINE
 470             size_type index1 () const {
 471                 return it1_;
 472             }
 473             BOOST_UBLAS_INLINE
 474             size_type index2 () const {
 475                 return it2_;
 476             }
 477
 478             // Assignment
 479             BOOST_UBLAS_INLINE
 480             const_iterator1 &operator = (const const_iterator1 &it) {
 481                 container_const_reference<self_type>::assign (&it ());
 482                 it1_ = it.it1_;
 483                 it2_ = it.it2_;
 484                 return *this;
 485             }
 486
 487             // Comparison
 488             BOOST_UBLAS_INLINE
 489             bool operator == (const const_iterator1 &it) const {
 490                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
 491                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
 492                 return it1_ == it.it1_;
 493             }
 494             BOOST_UBLAS_INLINE
 495             bool operator < (const const_iterator1 &it) const {
 496                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
 497                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
 498                 return it1_ < it.it1_;
 499             }
 500
 501         private:
 502             size_type it1_;
 503             size_type it2_;
 504         };
 505 #endif
 506
 507         BOOST_UBLAS_INLINE
 508         const_iterator1 begin1 () const {
 509             return find1 (0, 0, 0);
 510         }
 511         BOOST_UBLAS_INLINE
 512         const_iterator1 cbegin1 () const {
 513             return begin1 ();
 514         }
 515         BOOST_UBLAS_INLINE
 516         const_iterator1 end1 () const {
 517             return find1 (0, size1_, 0);
 518         }
 519         BOOST_UBLAS_INLINE
 520         const_iterator1 cend1 () const {
 521             return end1 ();
 522         }
 523
 524 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
 525         class iterator1:
 526             public container_reference<triangular_matrix>,
 527             public random_access_iterator_base<packed_random_access_iterator_tag,
 528                                                iterator1, value_type> {
 529         public:
 530             typedef typename triangular_matrix::value_type value_type;
 531             typedef typename triangular_matrix::difference_type difference_type;
 532             typedef typename triangular_matrix::reference reference;
 533             typedef typename triangular_matrix::pointer pointer;
 534
 535             typedef iterator2 dual_iterator_type;
 536             typedef reverse_iterator2 dual_reverse_iterator_type;
 537
 538             // Construction and destruction
 539             BOOST_UBLAS_INLINE
 540             iterator1 ():
 541                 container_reference<self_type> (), it1_ (), it2_ () {}
 542             BOOST_UBLAS_INLINE
 543             iterator1 (self_type &m, size_type it1, size_type it2):
 544                 container_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
 545
 546             // Arithmetic
 547             BOOST_UBLAS_INLINE
 548             iterator1 &operator ++ () {
 549                 ++ it1_;
 550                 return *this;
 551             }
 552             BOOST_UBLAS_INLINE
 553             iterator1 &operator -- () {
 554                 -- it1_;
 555                 return *this;
 556             }
 557             BOOST_UBLAS_INLINE
 558             iterator1 &operator += (difference_type n) {
 559                 it1_ += n;
 560                 return *this;
 561             }
 562             BOOST_UBLAS_INLINE
 563             iterator1 &operator -= (difference_type n) {
 564                 it1_ -= n;
 565                 return *this;
 566             }
 567             BOOST_UBLAS_INLINE
 568             difference_type operator - (const iterator1 &it) const {
 569                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
 570                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
 571                 return it1_ - it.it1_;
 572             }
 573
 574             // Dereference
 575             BOOST_UBLAS_INLINE
 576             reference operator * () const {
 577                 return (*this) () (it1_, it2_);
 578             }
 579             BOOST_UBLAS_INLINE
 580             reference operator [] (difference_type n) const {
 581                 return *(*this + n);
 582             }
 583
 584 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
 585             BOOST_UBLAS_INLINE
 586 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 587             typename self_type::
 588 #endif
 589             iterator2 begin () const {
 590                 return (*this) ().find2 (1, it1_, 0);
 591             }
 592             BOOST_UBLAS_INLINE
 593 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 594             typename self_type::
 595 #endif
 596             iterator2 end () const {
 597                 return (*this) ().find2 (1, it1_, (*this) ().size2 ());
 598             }
 599             BOOST_UBLAS_INLINE
 600 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 601             typename self_type::
 602 #endif
 603             reverse_iterator2 rbegin () const {
 604                 return reverse_iterator2 (end ());
 605             }
 606             BOOST_UBLAS_INLINE
 607 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 608             typename self_type::
 609 #endif
 610             reverse_iterator2 rend () const {
 611                 return reverse_iterator2 (begin ());
 612             }
 613 #endif
 614
 615             // Indices
 616             BOOST_UBLAS_INLINE
 617             size_type index1 () const {
 618                 return it1_;
 619             }
 620             BOOST_UBLAS_INLINE
 621             size_type index2 () const {
 622                 return it2_;
 623             }
 624
 625             // Assignment
 626             BOOST_UBLAS_INLINE
 627             iterator1 &operator = (const iterator1 &it) {
 628                 container_reference<self_type>::assign (&it ());
 629                 it1_ = it.it1_;
 630                 it2_ = it.it2_;
 631                 return *this;
 632             }
 633
 634             // Comparison
 635             BOOST_UBLAS_INLINE
 636             bool operator == (const iterator1 &it) const {
 637                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
 638                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
 639                 return it1_ == it.it1_;
 640             }
 641             BOOST_UBLAS_INLINE
 642             bool operator < (const iterator1 &it) const {
 643                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
 644                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
 645                 return it1_ < it.it1_;
 646             }
 647
 648         private:
 649             size_type it1_;
 650             size_type it2_;
 651
 652             friend class const_iterator1;
 653         };
 654 #endif
 655
 656         BOOST_UBLAS_INLINE
 657         iterator1 begin1 () {
 658             return find1 (0, 0, 0);
 659         }
 660         BOOST_UBLAS_INLINE
 661         iterator1 end1 () {
 662             return find1 (0, size1_, 0);
 663         }
 664
 665 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
 666         class const_iterator2:
 667             public container_const_reference<triangular_matrix>,
 668             public random_access_iterator_base<packed_random_access_iterator_tag,
 669                                                const_iterator2, value_type> {
 670         public:
 671             typedef typename triangular_matrix::value_type value_type;
 672             typedef typename triangular_matrix::difference_type difference_type;
 673             typedef typename triangular_matrix::const_reference reference;
 674             typedef const typename triangular_matrix::pointer pointer;
 675
 676             typedef const_iterator1 dual_iterator_type;
 677             typedef const_reverse_iterator1 dual_reverse_iterator_type;
 678
 679             // Construction and destruction
 680             BOOST_UBLAS_INLINE
 681             const_iterator2 ():
 682                 container_const_reference<self_type> (), it1_ (), it2_ () {}
 683             BOOST_UBLAS_INLINE
 684             const_iterator2 (const self_type &m, size_type it1, size_type it2):
 685                 container_const_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
 686             BOOST_UBLAS_INLINE
 687             const_iterator2 (const iterator2 &it):
 688                 container_const_reference<self_type> (it ()), it1_ (it.it1_), it2_ (it.it2_) {}
 689
 690             // Arithmetic
 691             BOOST_UBLAS_INLINE
 692             const_iterator2 &operator ++ () {
 693                 ++ it2_;
 694                 return *this;
 695             }
 696             BOOST_UBLAS_INLINE
 697             const_iterator2 &operator -- () {
 698                 -- it2_;
 699                 return *this;
 700             }
 701             BOOST_UBLAS_INLINE
 702             const_iterator2 &operator += (difference_type n) {
 703                 it2_ += n;
 704                 return *this;
 705             }
 706             BOOST_UBLAS_INLINE
 707             const_iterator2 &operator -= (difference_type n) {
 708                 it2_ -= n;
 709                 return *this;
 710             }
 711             BOOST_UBLAS_INLINE
 712             difference_type operator - (const const_iterator2 &it) const {
 713                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
 714                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
 715                 return it2_ - it.it2_;
 716             }
 717
 718             // Dereference
 719             BOOST_UBLAS_INLINE
 720             const_reference operator * () const {
 721                 return (*this) () (it1_, it2_);
 722             }
 723             BOOST_UBLAS_INLINE
 724             const_reference operator [] (difference_type n) const {
 725                 return *(*this + n);
 726             }
 727
 728 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
 729             BOOST_UBLAS_INLINE
 730 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 731             typename self_type::
 732 #endif
 733             const_iterator1 begin () const {
 734                 return (*this) ().find1 (1, 0, it2_);
 735             }
 736             BOOST_UBLAS_INLINE
 737 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 738             typename self_type::
 739 #endif
 740             const_iterator1 cbegin () const {
 741                 return begin ();
 742             }
 743             BOOST_UBLAS_INLINE
 744 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 745             typename self_type::
 746 #endif
 747             const_iterator1 end () const {
 748                 return (*this) ().find1 (1, (*this) ().size1 (), it2_);
 749             }
 750             BOOST_UBLAS_INLINE
 751 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 752             typename self_type::
 753 #endif
 754             const_iterator1 cend () const {
 755                 return end ();
 756             }
 757             BOOST_UBLAS_INLINE
 758 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 759             typename self_type::
 760 #endif
 761             const_reverse_iterator1 rbegin () const {
 762                 return const_reverse_iterator1 (end ());
 763             }
 764             BOOST_UBLAS_INLINE
 765 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 766             typename self_type::
 767 #endif
 768             const_reverse_iterator1 crbegin () const {
 769                 return rbegin ();
 770             }
 771
 772             BOOST_UBLAS_INLINE
 773 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 774             typename self_type::
 775 #endif
 776             const_reverse_iterator1 rend () const {
 777                 return const_reverse_iterator1 (begin ());
 778             }
 779             BOOST_UBLAS_INLINE
 780 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 781             typename self_type::
 782 #endif
 783             const_reverse_iterator1 crend () const {
 784                 return rend ();
 785             }
 786 #endif
 787
 788             // Indices
 789             BOOST_UBLAS_INLINE
 790             size_type index1 () const {
 791                 return it1_;
 792             }
 793             BOOST_UBLAS_INLINE
 794             size_type index2 () const {
 795                 return it2_;
 796             }
 797
 798             // Assignment
 799             BOOST_UBLAS_INLINE
 800             const_iterator2 &operator = (const const_iterator2 &it) {
 801                 container_const_reference<self_type>::assign (&it ());
 802                 it1_ = it.it1_;
 803                 it2_ = it.it2_;
 804                 return *this;
 805             }
 806
 807             // Comparison
 808             BOOST_UBLAS_INLINE
 809             bool operator == (const const_iterator2 &it) const {
 810                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
 811                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
 812                 return it2_ == it.it2_;
 813             }
 814             BOOST_UBLAS_INLINE
 815             bool operator < (const const_iterator2 &it) const {
 816                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
 817                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
 818                 return it2_ < it.it2_;
 819             }
 820
 821         private:
 822             size_type it1_;
 823             size_type it2_;
 824         };
 825 #endif
 826
 827         BOOST_UBLAS_INLINE
 828         const_iterator2 begin2 () const {
 829             return find2 (0, 0, 0);
 830         }
 831         BOOST_UBLAS_INLINE
 832         const_iterator2 cbegin2 () const {
 833             return begin2 ();
 834         }
 835         BOOST_UBLAS_INLINE
 836         const_iterator2 end2 () const {
 837             return find2 (0, 0, size2_);
 838         }
 839         BOOST_UBLAS_INLINE
 840         const_iterator2 cend2 () const {
 841             return end2 ();
 842         }
 843
 844 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
 845         class iterator2:
 846             public container_reference<triangular_matrix>,
 847             public random_access_iterator_base<packed_random_access_iterator_tag,
 848                                                iterator2, value_type> {
 849         public:
 850             typedef typename triangular_matrix::value_type value_type;
 851             typedef typename triangular_matrix::difference_type difference_type;
 852             typedef typename triangular_matrix::reference reference;
 853             typedef typename triangular_matrix::pointer pointer;
 854
 855             typedef iterator1 dual_iterator_type;
 856             typedef reverse_iterator1 dual_reverse_iterator_type;
 857
 858             // Construction and destruction
 859             BOOST_UBLAS_INLINE
 860             iterator2 ():
 861                 container_reference<self_type> (), it1_ (), it2_ () {}
 862             BOOST_UBLAS_INLINE
 863             iterator2 (self_type &m, size_type it1, size_type it2):
 864                 container_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
 865
 866             // Arithmetic
 867             BOOST_UBLAS_INLINE
 868             iterator2 &operator ++ () {
 869                 ++ it2_;
 870                 return *this;
 871             }
 872             BOOST_UBLAS_INLINE
 873             iterator2 &operator -- () {
 874                 -- it2_;
 875                 return *this;
 876             }
 877             BOOST_UBLAS_INLINE
 878             iterator2 &operator += (difference_type n) {
 879                 it2_ += n;
 880                 return *this;
 881             }
 882             BOOST_UBLAS_INLINE
 883             iterator2 &operator -= (difference_type n) {
 884                 it2_ -= n;
 885                 return *this;
 886             }
 887             BOOST_UBLAS_INLINE
 888             difference_type operator - (const iterator2 &it) const {
 889                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
 890                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
 891                 return it2_ - it.it2_;
 892             }
 893
 894             // Dereference
 895             BOOST_UBLAS_INLINE
 896             reference operator * () const {
 897                 return (*this) () (it1_, it2_);
 898             }
 899             BOOST_UBLAS_INLINE
 900             reference operator [] (difference_type n) const {
 901                 return *(*this + n);
 902             }
 903
 904 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
 905             BOOST_UBLAS_INLINE
 906 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 907             typename self_type::
 908 #endif
 909             iterator1 begin () const {
 910                 return (*this) ().find1 (1, 0, it2_);
 911             }
 912             BOOST_UBLAS_INLINE
 913 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 914             typename self_type::
 915 #endif
 916             iterator1 end () const {
 917                 return (*this) ().find1 (1, (*this) ().size1 (), it2_);
 918             }
 919             BOOST_UBLAS_INLINE
 920 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 921             typename self_type::
 922 #endif
 923             reverse_iterator1 rbegin () const {
 924                 return reverse_iterator1 (end ());
 925             }
 926             BOOST_UBLAS_INLINE
 927 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
 928             typename self_type::
 929 #endif
 930             reverse_iterator1 rend () const {
 931                 return reverse_iterator1 (begin ());
 932             }
 933 #endif
 934
 935             // Indices
 936             BOOST_UBLAS_INLINE
 937             size_type index1 () const {
 938                 return it1_;
 939             }
 940             BOOST_UBLAS_INLINE
 941             size_type index2 () const {
 942                 return it2_;
 943             }
 944
 945             // Assignment
 946             BOOST_UBLAS_INLINE
 947             iterator2 &operator = (const iterator2 &it) {
 948                 container_reference<self_type>::assign (&it ());
 949                 it1_ = it.it1_;
 950                 it2_ = it.it2_;
 951                 return *this;
 952             }
 953
 954             // Comparison
 955             BOOST_UBLAS_INLINE
 956             bool operator == (const iterator2 &it) const {
 957                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
 958                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
 959                 return it2_ == it.it2_;
 960             }
 961             BOOST_UBLAS_INLINE
 962             bool operator < (const iterator2 &it) const {
 963                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
 964                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
 965                 return it2_ < it.it2_;
 966             }
 967
 968         private:
 969             size_type it1_;
 970             size_type it2_;
 971
 972             friend class const_iterator2;
 973         };
 974 #endif
 975
 976         BOOST_UBLAS_INLINE
 977         iterator2 begin2 () {
 978             return find2 (0, 0, 0);
 979         }
 980         BOOST_UBLAS_INLINE
 981         iterator2 end2 () {
 982             return find2 (0, 0, size2_);
 983         }
 984
 985         // Reverse iterators
 986
 987         BOOST_UBLAS_INLINE
 988         const_reverse_iterator1 rbegin1 () const {
 989             return const_reverse_iterator1 (end1 ());
 990         }
 991         BOOST_UBLAS_INLINE
 992         const_reverse_iterator1 crbegin1 () const {
 993             return rbegin1 ();
 994         }
 995         BOOST_UBLAS_INLINE
 996         const_reverse_iterator1 rend1 () const {
 997             return const_reverse_iterator1 (begin1 ());
 998         }
 999         BOOST_UBLAS_INLINE
1000         const_reverse_iterator1 crend1 () const {
1001             return rend1 ();
1002         }
1003
1004         BOOST_UBLAS_INLINE
1005         reverse_iterator1 rbegin1 () {
1006             return reverse_iterator1 (end1 ());
1007         }
1008         BOOST_UBLAS_INLINE
1009         reverse_iterator1 rend1 () {
1010             return reverse_iterator1 (begin1 ());
1011         }
1012
1013         BOOST_UBLAS_INLINE
1014         const_reverse_iterator2 rbegin2 () const {
1015             return const_reverse_iterator2 (end2 ());
1016         }
1017         BOOST_UBLAS_INLINE
1018         const_reverse_iterator2 crbegin2 () const {
1019             return rbegin2 ();
1020         }
1021         BOOST_UBLAS_INLINE
1022         const_reverse_iterator2 rend2 () const {
1023             return const_reverse_iterator2 (begin2 ());
1024         }
1025         BOOST_UBLAS_INLINE
1026         const_reverse_iterator2 crend2 () const {
1027             return rend2 ();
1028         }
1029
1030         BOOST_UBLAS_INLINE
1031         reverse_iterator2 rbegin2 () {
1032             return reverse_iterator2 (end2 ());
1033         }
1034         BOOST_UBLAS_INLINE
1035         reverse_iterator2 rend2 () {
1036             return reverse_iterator2 (begin2 ());
1037         }
1038
1039     private:
1040         size_type size1_;
1041         size_type size2_;
1042         array_type data_;
1043         static const value_type zero_;
1044         static const value_type one_;
1045     };
1046
1047     template<class T, class TRI, class L, class A>
1048     const typename triangular_matrix<T, TRI, L, A>::value_type triangular_matrix<T, TRI, L, A>::zero_ = value_type/*zero*/();
1049     template<class T, class TRI, class L, class A>
1050     const typename triangular_matrix<T, TRI, L, A>::value_type triangular_matrix<T, TRI, L, A>::one_ (1);
1051
1052
1053     // Triangular matrix adaptor class
1054     template<class M, class TRI>
1055     class triangular_adaptor:
1056         public matrix_expression<triangular_adaptor<M, TRI> > {
1057
1058         typedef triangular_adaptor<M, TRI> self_type;
1059
1060     public:
1061 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
1062         using matrix_expression<self_type>::operator ();
1063 #endif
1064         typedef const M const_matrix_type;
1065         typedef M matrix_type;
1066         typedef TRI triangular_type;
1067         typedef typename M::size_type size_type;
1068         typedef typename M::difference_type difference_type;
1069         typedef typename M::value_type value_type;
1070         typedef typename M::const_reference const_reference;
1071         typedef typename boost::mpl::if_<boost::is_const<M>,
1072                                           typename M::const_reference,
1073                                           typename M::reference>::type reference;
1074         typedef typename boost::mpl::if_<boost::is_const<M>,
1075                                           typename M::const_closure_type,
1076                                           typename M::closure_type>::type matrix_closure_type;
1077         typedef const self_type const_closure_type;
1078         typedef self_type closure_type;
1079         // Replaced by _temporary_traits to avoid type requirements on M
1080         //typedef typename M::vector_temporary_type vector_temporary_type;
1081         //typedef typename M::matrix_temporary_type matrix_temporary_type;
1082         typedef typename storage_restrict_traits<typename M::storage_category,
1083                                                  packed_proxy_tag>::storage_category storage_category;
1084         typedef typename M::orientation_category orientation_category;
1085
1086         // Construction and destruction
1087         BOOST_UBLAS_INLINE
1088         triangular_adaptor (matrix_type &data):
1089             matrix_expression<self_type> (),
1090             data_ (data) {}
1091         BOOST_UBLAS_INLINE
1092         triangular_adaptor (const triangular_adaptor &m):
1093             matrix_expression<self_type> (),
1094             data_ (m.data_) {}
1095
1096         // Accessors
1097         BOOST_UBLAS_INLINE
1098         size_type size1 () const {
1099             return data_.size1 ();
1100         }
1101         BOOST_UBLAS_INLINE
1102         size_type size2 () const {
1103             return data_.size2 ();
1104         }
1105
1106         // Storage accessors
1107         BOOST_UBLAS_INLINE
1108         const matrix_closure_type &data () const {
1109             return data_;
1110         }
1111         BOOST_UBLAS_INLINE
1112         matrix_closure_type &data () {
1113             return data_;
1114         }
1115
1116         // Element access
1117 #ifndef BOOST_UBLAS_PROXY_CONST_MEMBER
1118         BOOST_UBLAS_INLINE
1119         const_reference operator () (size_type i, size_type j) const {
1120             BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
1121             BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
1122             if (triangular_type::other (i, j))
1123                 return data () (i, j);
1124             else if (triangular_type::one (i, j))
1125                 return one_;
1126             else
1127                 return zero_;
1128         }
1129         BOOST_UBLAS_INLINE
1130         reference operator () (size_type i, size_type j) {
1131             BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
1132             BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
1133             if (!triangular_type::other (i, j)) {
1134                 bad_index ().raise ();
1135                 // NEVER reached
1136             }
1137             return data () (i, j);
1138         }
1139 #else
1140         BOOST_UBLAS_INLINE
1141         reference operator () (size_type i, size_type j) const {
1142             BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
1143             BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
1144             if (!triangular_type::other (i, j)) {
1145                 bad_index ().raise ();
1146                 // NEVER reached
1147             }
1148             return data () (i, j);
1149         }
1150 #endif
1151
1152         // Assignment
1153         BOOST_UBLAS_INLINE
1154         triangular_adaptor &operator = (const triangular_adaptor &m) {
1155             matrix_assign<scalar_assign> (*this, m);
1156             return *this;
1157         }
1158         BOOST_UBLAS_INLINE
1159         triangular_adaptor &assign_temporary (triangular_adaptor &m) {
1160             *this = m;
1161             return *this;
1162         }
1163         template<class AE>
1164         BOOST_UBLAS_INLINE
1165         triangular_adaptor &operator = (const matrix_expression<AE> &ae) {
1166             matrix_assign<scalar_assign> (*this, matrix<value_type> (ae));
1167             return *this;
1168         }
1169         template<class AE>
1170         BOOST_UBLAS_INLINE
1171         triangular_adaptor &assign (const matrix_expression<AE> &ae) {
1172             matrix_assign<scalar_assign> (*this, ae);
1173             return *this;
1174         }
1175         template<class AE>
1176         BOOST_UBLAS_INLINE
1177         triangular_adaptor& operator += (const matrix_expression<AE> &ae) {
1178             matrix_assign<scalar_assign> (*this, matrix<value_type> (*this + ae));
1179             return *this;
1180         }
1181         template<class AE>
1182         BOOST_UBLAS_INLINE
1183         triangular_adaptor &plus_assign (const matrix_expression<AE> &ae) {
1184             matrix_assign<scalar_plus_assign> (*this, ae);
1185             return *this;
1186         }
1187         template<class AE>
1188         BOOST_UBLAS_INLINE
1189         triangular_adaptor& operator -= (const matrix_expression<AE> &ae) {
1190             matrix_assign<scalar_assign> (*this, matrix<value_type> (*this - ae));
1191             return *this;
1192         }
1193         template<class AE>
1194         BOOST_UBLAS_INLINE
1195         triangular_adaptor &minus_assign (const matrix_expression<AE> &ae) {
1196             matrix_assign<scalar_minus_assign> (*this, ae);
1197             return *this;
1198         }
1199         template<class AT>
1200         BOOST_UBLAS_INLINE
1201         triangular_adaptor& operator *= (const AT &at) {
1202             matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
1203             return *this;
1204         }
1205         template<class AT>
1206         BOOST_UBLAS_INLINE
1207         triangular_adaptor& operator /= (const AT &at) {
1208             matrix_assign_scalar<scalar_divides_assign> (*this, at);
1209             return *this;
1210         }
1211
1212         // Closure comparison
1213         BOOST_UBLAS_INLINE
1214         bool same_closure (const triangular_adaptor &ta) const {
1215             return (*this).data ().same_closure (ta.data ());
1216         }
1217
1218         // Swapping
1219         BOOST_UBLAS_INLINE
1220         void swap (triangular_adaptor &m) {
1221             if (this != &m)
1222                 matrix_swap<scalar_swap> (*this, m);
1223         }
1224         BOOST_UBLAS_INLINE
1225         friend void swap (triangular_adaptor &m1, triangular_adaptor &m2) {
1226             m1.swap (m2);
1227         }
1228
1229         // Iterator types
1230    private:
1231         typedef typename M::const_iterator1 const_subiterator1_type;
1232         typedef typename boost::mpl::if_<boost::is_const<M>,
1233                                           typename M::const_iterator1,
1234                                           typename M::iterator1>::type subiterator1_type;
1235         typedef typename M::const_iterator2 const_subiterator2_type;
1236         typedef typename boost::mpl::if_<boost::is_const<M>,
1237                                           typename M::const_iterator2,
1238                                           typename M::iterator2>::type subiterator2_type;
1239
1240     public:
1241 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1242         typedef indexed_iterator1<self_type, packed_random_access_iterator_tag> iterator1;
1243         typedef indexed_iterator2<self_type, packed_random_access_iterator_tag> iterator2;
1244         typedef indexed_const_iterator1<self_type, packed_random_access_iterator_tag> const_iterator1;
1245         typedef indexed_const_iterator2<self_type, packed_random_access_iterator_tag> const_iterator2;
1246 #else
1247         class const_iterator1;
1248         class iterator1;
1249         class const_iterator2;
1250         class iterator2;
1251 #endif
1252         typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
1253         typedef reverse_iterator_base1<iterator1> reverse_iterator1;
1254         typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
1255         typedef reverse_iterator_base2<iterator2> reverse_iterator2;
1256
1257         // Element lookup
1258         BOOST_UBLAS_INLINE
1259         const_iterator1 find1 (int rank, size_type i, size_type j) const {
1260             if (rank == 1)
1261                 i = triangular_type::restrict1 (i, j, size1(), size2());
1262             if (rank == 0)
1263                 i = triangular_type::global_restrict1 (i, size1(), j, size2());
1264             return const_iterator1 (*this, data ().find1 (rank, i, j));
1265         }
1266         BOOST_UBLAS_INLINE
1267         iterator1 find1 (int rank, size_type i, size_type j) {
1268             if (rank == 1)
1269                 i = triangular_type::mutable_restrict1 (i, j, size1(), size2());
1270             if (rank == 0)
1271                 i = triangular_type::global_mutable_restrict1 (i, size1(), j, size2());
1272             return iterator1 (*this, data ().find1 (rank, i, j));
1273         }
1274         BOOST_UBLAS_INLINE
1275         const_iterator2 find2 (int rank, size_type i, size_type j) const {
1276             if (rank == 1)
1277                 j = triangular_type::restrict2 (i, j, size1(), size2());
1278             if (rank == 0)
1279                 j = triangular_type::global_restrict2 (i, size1(), j, size2());
1280             return const_iterator2 (*this, data ().find2 (rank, i, j));
1281         }
1282         BOOST_UBLAS_INLINE
1283         iterator2 find2 (int rank, size_type i, size_type j) {
1284             if (rank == 1)
1285                 j = triangular_type::mutable_restrict2 (i, j, size1(), size2());
1286             if (rank == 0)
1287                 j = triangular_type::global_mutable_restrict2 (i, size1(), j, size2());
1288             return iterator2 (*this, data ().find2 (rank, i, j));
1289         }
1290
1291         // Iterators simply are indices.
1292
1293 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1294         class const_iterator1:
1295             public container_const_reference<triangular_adaptor>,
1296             public random_access_iterator_base<typename iterator_restrict_traits<
1297                                                    typename const_subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,
1298                                                const_iterator1, value_type> {
1299         public:
1300             typedef typename const_subiterator1_type::value_type value_type;
1301             typedef typename const_subiterator1_type::difference_type difference_type;
1302             typedef typename const_subiterator1_type::reference reference;
1303             typedef typename const_subiterator1_type::pointer pointer;
1304
1305             typedef const_iterator2 dual_iterator_type;
1306             typedef const_reverse_iterator2 dual_reverse_iterator_type;
1307
1308             // Construction and destruction
1309             BOOST_UBLAS_INLINE
1310             const_iterator1 ():
1311                 container_const_reference<self_type> (), it1_ () {}
1312             BOOST_UBLAS_INLINE
1313             const_iterator1 (const self_type &m, const const_subiterator1_type &it1):
1314                 container_const_reference<self_type> (m), it1_ (it1) {}
1315             BOOST_UBLAS_INLINE
1316             const_iterator1 (const iterator1 &it):
1317                 container_const_reference<self_type> (it ()), it1_ (it.it1_) {}
1318
1319             // Arithmetic
1320             BOOST_UBLAS_INLINE
1321             const_iterator1 &operator ++ () {
1322                 ++ it1_;
1323                 return *this;
1324             }
1325             BOOST_UBLAS_INLINE
1326             const_iterator1 &operator -- () {
1327                 -- it1_;
1328                 return *this;
1329             }
1330             BOOST_UBLAS_INLINE
1331             const_iterator1 &operator += (difference_type n) {
1332                 it1_ += n;
1333                 return *this;
1334             }
1335             BOOST_UBLAS_INLINE
1336             const_iterator1 &operator -= (difference_type n) {
1337                 it1_ -= n;
1338                 return *this;
1339             }
1340             BOOST_UBLAS_INLINE
1341             difference_type operator - (const const_iterator1 &it) const {
1342                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1343                 return it1_ - it.it1_;
1344             }
1345
1346             // Dereference
1347             BOOST_UBLAS_INLINE
1348             const_reference operator * () const {
1349                 size_type i = index1 ();
1350                 size_type j = index2 ();
1351                 BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());
1352                 BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());
1353                 if (triangular_type::other (i, j))
1354                     return *it1_;
1355                 else
1356                     return (*this) () (i, j);
1357             }
1358             BOOST_UBLAS_INLINE
1359             const_reference operator [] (difference_type n) const {
1360                 return *(*this + n);
1361             }
1362
1363 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1364             BOOST_UBLAS_INLINE
1365 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1366             typename self_type::
1367 #endif
1368             const_iterator2 begin () const {
1369                 return (*this) ().find2 (1, index1 (), 0);
1370             }
1371             BOOST_UBLAS_INLINE
1372 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1373             typename self_type::
1374 #endif
1375             const_iterator2 cbegin () const {
1376                 return begin ();
1377             }
1378             BOOST_UBLAS_INLINE
1379 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1380             typename self_type::
1381 #endif
1382             const_iterator2 end () const {
1383                 return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
1384             }
1385             BOOST_UBLAS_INLINE
1386 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1387             typename self_type::
1388 #endif
1389             const_iterator2 cend () const {
1390                 return end ();
1391             }
1392
1393             BOOST_UBLAS_INLINE
1394 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1395             typename self_type::
1396 #endif
1397             const_reverse_iterator2 rbegin () const {
1398                 return const_reverse_iterator2 (end ());
1399             }
1400             BOOST_UBLAS_INLINE
1401 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1402             typename self_type::
1403 #endif
1404             const_reverse_iterator2 crbegin () const {
1405                 return rbegin ();
1406             }
1407             BOOST_UBLAS_INLINE
1408 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1409             typename self_type::
1410 #endif
1411             const_reverse_iterator2 rend () const {
1412                 return const_reverse_iterator2 (begin ());
1413             }
1414             BOOST_UBLAS_INLINE
1415 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1416             typename self_type::
1417 #endif
1418             const_reverse_iterator2 crend () const {
1419                 return rend ();
1420             }
1421 #endif
1422
1423             // Indices
1424             BOOST_UBLAS_INLINE
1425             size_type index1 () const {
1426                 return it1_.index1 ();
1427             }
1428             BOOST_UBLAS_INLINE
1429             size_type index2 () const {
1430                 return it1_.index2 ();
1431             }
1432
1433             // Assignment
1434             BOOST_UBLAS_INLINE
1435             const_iterator1 &operator = (const const_iterator1 &it) {
1436                 container_const_reference<self_type>::assign (&it ());
1437                 it1_ = it.it1_;
1438                 return *this;
1439             }
1440
1441             // Comparison
1442             BOOST_UBLAS_INLINE
1443             bool operator == (const const_iterator1 &it) const {
1444                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1445                 return it1_ == it.it1_;
1446             }
1447             BOOST_UBLAS_INLINE
1448             bool operator < (const const_iterator1 &it) const {
1449                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1450                 return it1_ < it.it1_;
1451             }
1452
1453         private:
1454             const_subiterator1_type it1_;
1455         };
1456 #endif
1457
1458         BOOST_UBLAS_INLINE
1459         const_iterator1 begin1 () const {
1460             return find1 (0, 0, 0);
1461         }
1462         BOOST_UBLAS_INLINE
1463         const_iterator1 cbegin1 () const {
1464             return begin1 ();
1465         }
1466         BOOST_UBLAS_INLINE
1467         const_iterator1 end1 () const {
1468             return find1 (0, size1 (), 0);
1469         }
1470         BOOST_UBLAS_INLINE
1471         const_iterator1 cend1 () const {
1472             return end1 ();
1473         }
1474
1475 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1476         class iterator1:
1477             public container_reference<triangular_adaptor>,
1478             public random_access_iterator_base<typename iterator_restrict_traits<
1479                                                    typename subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,
1480                                                iterator1, value_type> {
1481         public:
1482             typedef typename subiterator1_type::value_type value_type;
1483             typedef typename subiterator1_type::difference_type difference_type;
1484             typedef typename subiterator1_type::reference reference;
1485             typedef typename subiterator1_type::pointer pointer;
1486
1487             typedef iterator2 dual_iterator_type;
1488             typedef reverse_iterator2 dual_reverse_iterator_type;
1489
1490             // Construction and destruction
1491             BOOST_UBLAS_INLINE
1492             iterator1 ():
1493                 container_reference<self_type> (), it1_ () {}
1494             BOOST_UBLAS_INLINE
1495             iterator1 (self_type &m, const subiterator1_type &it1):
1496                 container_reference<self_type> (m), it1_ (it1) {}
1497
1498             // Arithmetic
1499             BOOST_UBLAS_INLINE
1500             iterator1 &operator ++ () {
1501                 ++ it1_;
1502                 return *this;
1503             }
1504             BOOST_UBLAS_INLINE
1505             iterator1 &operator -- () {
1506                 -- it1_;
1507                 return *this;
1508             }
1509             BOOST_UBLAS_INLINE
1510             iterator1 &operator += (difference_type n) {
1511                 it1_ += n;
1512                 return *this;
1513             }
1514             BOOST_UBLAS_INLINE
1515             iterator1 &operator -= (difference_type n) {
1516                 it1_ -= n;
1517                 return *this;
1518             }
1519             BOOST_UBLAS_INLINE
1520             difference_type operator - (const iterator1 &it) const {
1521                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1522                 return it1_ - it.it1_;
1523             }
1524
1525             // Dereference
1526             BOOST_UBLAS_INLINE
1527             reference operator * () const {
1528                 size_type i = index1 ();
1529                 size_type j = index2 ();
1530                 BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());
1531                 BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());
1532                 if (triangular_type::other (i, j))
1533                     return *it1_;
1534                 else
1535                     return (*this) () (i, j);
1536             }
1537             BOOST_UBLAS_INLINE
1538             reference operator [] (difference_type n) const {
1539                 return *(*this + n);
1540             }
1541
1542 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1543             BOOST_UBLAS_INLINE
1544 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1545             typename self_type::
1546 #endif
1547             iterator2 begin () const {
1548                 return (*this) ().find2 (1, index1 (), 0);
1549             }
1550             BOOST_UBLAS_INLINE
1551 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1552             typename self_type::
1553 #endif
1554             iterator2 end () const {
1555                 return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
1556             }
1557             BOOST_UBLAS_INLINE
1558 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1559             typename self_type::
1560 #endif
1561             reverse_iterator2 rbegin () const {
1562                 return reverse_iterator2 (end ());
1563             }
1564             BOOST_UBLAS_INLINE
1565 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1566             typename self_type::
1567 #endif
1568             reverse_iterator2 rend () const {
1569                 return reverse_iterator2 (begin ());
1570             }
1571 #endif
1572
1573             // Indices
1574             BOOST_UBLAS_INLINE
1575             size_type index1 () const {
1576                 return it1_.index1 ();
1577             }
1578             BOOST_UBLAS_INLINE
1579             size_type index2 () const {
1580                 return it1_.index2 ();
1581             }
1582
1583             // Assignment
1584             BOOST_UBLAS_INLINE
1585             iterator1 &operator = (const iterator1 &it) {
1586                 container_reference<self_type>::assign (&it ());
1587                 it1_ = it.it1_;
1588                 return *this;
1589             }
1590
1591             // Comparison
1592             BOOST_UBLAS_INLINE
1593             bool operator == (const iterator1 &it) const {
1594                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1595                 return it1_ == it.it1_;
1596             }
1597             BOOST_UBLAS_INLINE
1598             bool operator < (const iterator1 &it) const {
1599                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1600                 return it1_ < it.it1_;
1601             }
1602
1603         private:
1604             subiterator1_type it1_;
1605
1606             friend class const_iterator1;
1607         };
1608 #endif
1609
1610         BOOST_UBLAS_INLINE
1611         iterator1 begin1 () {
1612             return find1 (0, 0, 0);
1613         }
1614         BOOST_UBLAS_INLINE
1615         iterator1 end1 () {
1616             return find1 (0, size1 (), 0);
1617         }
1618
1619 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1620         class const_iterator2:
1621             public container_const_reference<triangular_adaptor>,
1622             public random_access_iterator_base<typename iterator_restrict_traits<
1623                                                    typename const_subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,
1624                                                const_iterator2, value_type> {
1625         public:
1626             typedef typename const_subiterator2_type::value_type value_type;
1627             typedef typename const_subiterator2_type::difference_type difference_type;
1628             typedef typename const_subiterator2_type::reference reference;
1629             typedef typename const_subiterator2_type::pointer pointer;
1630
1631             typedef const_iterator1 dual_iterator_type;
1632             typedef const_reverse_iterator1 dual_reverse_iterator_type;
1633
1634             // Construction and destruction
1635             BOOST_UBLAS_INLINE
1636             const_iterator2 ():
1637                 container_const_reference<self_type> (), it2_ () {}
1638             BOOST_UBLAS_INLINE
1639             const_iterator2 (const self_type &m, const const_subiterator2_type &it2):
1640                 container_const_reference<self_type> (m), it2_ (it2) {}
1641             BOOST_UBLAS_INLINE
1642             const_iterator2 (const iterator2 &it):
1643                 container_const_reference<self_type> (it ()), it2_ (it.it2_) {}
1644
1645             // Arithmetic
1646             BOOST_UBLAS_INLINE
1647             const_iterator2 &operator ++ () {
1648                 ++ it2_;
1649                 return *this;
1650             }
1651             BOOST_UBLAS_INLINE
1652             const_iterator2 &operator -- () {
1653                 -- it2_;
1654                 return *this;
1655             }
1656             BOOST_UBLAS_INLINE
1657             const_iterator2 &operator += (difference_type n) {
1658                 it2_ += n;
1659                 return *this;
1660             }
1661             BOOST_UBLAS_INLINE
1662             const_iterator2 &operator -= (difference_type n) {
1663                 it2_ -= n;
1664                 return *this;
1665             }
1666             BOOST_UBLAS_INLINE
1667             difference_type operator - (const const_iterator2 &it) const {
1668                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1669                 return it2_ - it.it2_;
1670             }
1671
1672             // Dereference
1673             BOOST_UBLAS_INLINE
1674             const_reference operator * () const {
1675                 size_type i = index1 ();
1676                 size_type j = index2 ();
1677                 BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());
1678                 BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());
1679                 if (triangular_type::other (i, j))
1680                     return *it2_;
1681                 else
1682                     return (*this) () (i, j);
1683             }
1684             BOOST_UBLAS_INLINE
1685             const_reference operator [] (difference_type n) const {
1686                 return *(*this + n);
1687             }
1688
1689 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1690             BOOST_UBLAS_INLINE
1691 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1692             typename self_type::
1693 #endif
1694             const_iterator1 begin () const {
1695                 return (*this) ().find1 (1, 0, index2 ());
1696             }
1697             BOOST_UBLAS_INLINE
1698 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1699             typename self_type::
1700 #endif
1701             const_iterator1 cbegin () const {
1702                 return begin ();
1703             }
1704             BOOST_UBLAS_INLINE
1705 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1706             typename self_type::
1707 #endif
1708             const_iterator1 end () const {
1709                 return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
1710             }
1711             BOOST_UBLAS_INLINE
1712 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1713             typename self_type::
1714 #endif
1715             const_iterator1 cend () const {
1716                 return end ();
1717             }
1718             BOOST_UBLAS_INLINE
1719 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1720             typename self_type::
1721 #endif
1722             const_reverse_iterator1 rbegin () const {
1723                 return const_reverse_iterator1 (end ());
1724             }
1725             BOOST_UBLAS_INLINE
1726 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1727             typename self_type::
1728 #endif
1729             const_reverse_iterator1 crbegin () const {
1730                 return rbegin ();
1731             }
1732             BOOST_UBLAS_INLINE
1733 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1734             typename self_type::
1735 #endif
1736             const_reverse_iterator1 rend () const {
1737                 return const_reverse_iterator1 (begin ());
1738             }
1739             BOOST_UBLAS_INLINE
1740 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1741             typename self_type::
1742 #endif
1743             const_reverse_iterator1 crend () const {
1744                 return rend ();
1745             }
1746 #endif
1747
1748             // Indices
1749             BOOST_UBLAS_INLINE
1750             size_type index1 () const {
1751                 return it2_.index1 ();
1752             }
1753             BOOST_UBLAS_INLINE
1754             size_type index2 () const {
1755                 return it2_.index2 ();
1756             }
1757
1758             // Assignment
1759             BOOST_UBLAS_INLINE
1760             const_iterator2 &operator = (const const_iterator2 &it) {
1761                 container_const_reference<self_type>::assign (&it ());
1762                 it2_ = it.it2_;
1763                 return *this;
1764             }
1765
1766             // Comparison
1767             BOOST_UBLAS_INLINE
1768             bool operator == (const const_iterator2 &it) const {
1769                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1770                 return it2_ == it.it2_;
1771             }
1772             BOOST_UBLAS_INLINE
1773             bool operator < (const const_iterator2 &it) const {
1774                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1775                 return it2_ < it.it2_;
1776             }
1777
1778         private:
1779             const_subiterator2_type it2_;
1780         };
1781 #endif
1782
1783         BOOST_UBLAS_INLINE
1784         const_iterator2 begin2 () const {
1785             return find2 (0, 0, 0);
1786         }
1787         BOOST_UBLAS_INLINE
1788         const_iterator2 cbegin2 () const {
1789             return begin2 ();
1790         }
1791         BOOST_UBLAS_INLINE
1792         const_iterator2 end2 () const {
1793             return find2 (0, 0, size2 ());
1794         }
1795         BOOST_UBLAS_INLINE
1796         const_iterator2 cend2 () const {
1797             return end2 ();
1798         }
1799
1800 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1801         class iterator2:
1802             public container_reference<triangular_adaptor>,
1803             public random_access_iterator_base<typename iterator_restrict_traits<
1804                                                    typename subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,
1805                                                iterator2, value_type> {
1806         public:
1807             typedef typename subiterator2_type::value_type value_type;
1808             typedef typename subiterator2_type::difference_type difference_type;
1809             typedef typename subiterator2_type::reference reference;
1810             typedef typename subiterator2_type::pointer pointer;
1811
1812             typedef iterator1 dual_iterator_type;
1813             typedef reverse_iterator1 dual_reverse_iterator_type;
1814
1815             // Construction and destruction
1816             BOOST_UBLAS_INLINE
1817             iterator2 ():
1818                 container_reference<self_type> (), it2_ () {}
1819             BOOST_UBLAS_INLINE
1820             iterator2 (self_type &m, const subiterator2_type &it2):
1821                 container_reference<self_type> (m), it2_ (it2) {}
1822
1823             // Arithmetic
1824             BOOST_UBLAS_INLINE
1825             iterator2 &operator ++ () {
1826                 ++ it2_;
1827                 return *this;
1828             }
1829             BOOST_UBLAS_INLINE
1830             iterator2 &operator -- () {
1831                 -- it2_;
1832                 return *this;
1833             }
1834             BOOST_UBLAS_INLINE
1835             iterator2 &operator += (difference_type n) {
1836                 it2_ += n;
1837                 return *this;
1838             }
1839             BOOST_UBLAS_INLINE
1840             iterator2 &operator -= (difference_type n) {
1841                 it2_ -= n;
1842                 return *this;
1843             }
1844             BOOST_UBLAS_INLINE
1845             difference_type operator - (const iterator2 &it) const {
1846                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1847                 return it2_ - it.it2_;
1848             }
1849
1850             // Dereference
1851             BOOST_UBLAS_INLINE
1852             reference operator * () const {
1853                 size_type i = index1 ();
1854                 size_type j = index2 ();
1855                 BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());
1856                 BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());
1857                 if (triangular_type::other (i, j))
1858                     return *it2_;
1859                 else
1860                     return (*this) () (i, j);
1861             }
1862             BOOST_UBLAS_INLINE
1863             reference operator [] (difference_type n) const {
1864                 return *(*this + n);
1865             }
1866
1867 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1868             BOOST_UBLAS_INLINE
1869 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1870             typename self_type::
1871 #endif
1872             iterator1 begin () const {
1873                 return (*this) ().find1 (1, 0, index2 ());
1874             }
1875             BOOST_UBLAS_INLINE
1876 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1877             typename self_type::
1878 #endif
1879             iterator1 end () const {
1880                 return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
1881             }
1882             BOOST_UBLAS_INLINE
1883 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1884             typename self_type::
1885 #endif
1886             reverse_iterator1 rbegin () const {
1887                 return reverse_iterator1 (end ());
1888             }
1889             BOOST_UBLAS_INLINE
1890 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1891             typename self_type::
1892 #endif
1893             reverse_iterator1 rend () const {
1894                 return reverse_iterator1 (begin ());
1895             }
1896 #endif
1897
1898             // Indices
1899             BOOST_UBLAS_INLINE
1900             size_type index1 () const {
1901                 return it2_.index1 ();
1902             }
1903             BOOST_UBLAS_INLINE
1904             size_type index2 () const {
1905                 return it2_.index2 ();
1906             }
1907
1908             // Assignment
1909             BOOST_UBLAS_INLINE
1910             iterator2 &operator = (const iterator2 &it) {
1911                 container_reference<self_type>::assign (&it ());
1912                 it2_ = it.it2_;
1913                 return *this;
1914             }
1915
1916             // Comparison
1917             BOOST_UBLAS_INLINE
1918             bool operator == (const iterator2 &it) const {
1919                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1920                 return it2_ == it.it2_;
1921             }
1922             BOOST_UBLAS_INLINE
1923             bool operator < (const iterator2 &it) const {
1924                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1925                 return it2_ < it.it2_;
1926             }
1927
1928         private:
1929             subiterator2_type it2_;
1930
1931             friend class const_iterator2;
1932         };
1933 #endif
1934
1935         BOOST_UBLAS_INLINE
1936         iterator2 begin2 () {
1937             return find2 (0, 0, 0);
1938         }
1939         BOOST_UBLAS_INLINE
1940         iterator2 end2 () {
1941             return find2 (0, 0, size2 ());
1942         }
1943
1944         // Reverse iterators
1945
1946         BOOST_UBLAS_INLINE
1947         const_reverse_iterator1 rbegin1 () const {
1948             return const_reverse_iterator1 (end1 ());
1949         }
1950         BOOST_UBLAS_INLINE
1951         const_reverse_iterator1 crbegin1 () const {
1952             return rbegin1 ();
1953         }
1954         BOOST_UBLAS_INLINE
1955         const_reverse_iterator1 rend1 () const {
1956             return const_reverse_iterator1 (begin1 ());
1957         }
1958         BOOST_UBLAS_INLINE
1959         const_reverse_iterator1 crend1 () const {
1960             return rend1 ();
1961         }
1962
1963         BOOST_UBLAS_INLINE
1964         reverse_iterator1 rbegin1 () {
1965             return reverse_iterator1 (end1 ());
1966         }
1967         BOOST_UBLAS_INLINE
1968         reverse_iterator1 rend1 () {
1969             return reverse_iterator1 (begin1 ());
1970         }
1971
1972         BOOST_UBLAS_INLINE
1973         const_reverse_iterator2 rbegin2 () const {
1974             return const_reverse_iterator2 (end2 ());
1975         }
1976         BOOST_UBLAS_INLINE
1977         const_reverse_iterator2 crbegin2 () const {
1978             return rbegin2 ();
1979         }
1980         BOOST_UBLAS_INLINE
1981         const_reverse_iterator2 rend2 () const {
1982             return const_reverse_iterator2 (begin2 ());
1983         }
1984         BOOST_UBLAS_INLINE
1985         const_reverse_iterator2 crend2 () const {
1986             return rend2 ();
1987         }
1988
1989         BOOST_UBLAS_INLINE
1990         reverse_iterator2 rbegin2 () {
1991             return reverse_iterator2 (end2 ());
1992         }
1993         BOOST_UBLAS_INLINE
1994         reverse_iterator2 rend2 () {
1995             return reverse_iterator2 (begin2 ());
1996         }
1997
1998     private:
1999         matrix_closure_type data_;
2000         static const value_type zero_;
2001         static const value_type one_;
2002     };
2003
2004     template<class M, class TRI>
2005     const typename triangular_adaptor<M, TRI>::value_type triangular_adaptor<M, TRI>::zero_ = value_type/*zero*/();
2006     template<class M, class TRI>
2007     const typename triangular_adaptor<M, TRI>::value_type triangular_adaptor<M, TRI>::one_ (1);
2008
2009     template <class M, class TRI>
2010     struct vector_temporary_traits< triangular_adaptor<M, TRI> >
2011     : vector_temporary_traits< typename boost::remove_const<M>::type > {} ;
2012     template <class M, class TRI>
2013     struct vector_temporary_traits< const triangular_adaptor<M, TRI> >
2014     : vector_temporary_traits< typename boost::remove_const<M>::type > {} ;
2015
2016     template <class M, class TRI>
2017     struct matrix_temporary_traits< triangular_adaptor<M, TRI> >
2018     : matrix_temporary_traits< typename boost::remove_const<M>::type > {};
2019     template <class M, class TRI>
2020     struct matrix_temporary_traits< const triangular_adaptor<M, TRI> >
2021     : matrix_temporary_traits< typename boost::remove_const<M>::type > {};
2022
2023
2024     template<class E1, class E2>
2025     struct matrix_vector_solve_traits {
2026         typedef typename promote_traits<typename E1::value_type, typename E2::value_type>::promote_type promote_type;
2027         typedef vector<promote_type> result_type;
2028     };
2029
2030     // Operations:
2031     //  n * (n - 1) / 2 + n = n * (n + 1) / 2 multiplications,
2032     //  n * (n - 1) / 2 additions
2033
2034     // Dense (proxy) case
2035     template<class E1, class E2>
2036     BOOST_UBLAS_INLINE
2037     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
2038                         lower_tag, column_major_tag, dense_proxy_tag) {
2039         typedef typename E2::size_type size_type;
2040         typedef typename E2::value_type value_type;
2041
2042         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
2043         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
2044         size_type size = e2 ().size ();
2045         for (size_type n = 0; n < size; ++ n) {
2046 #ifndef BOOST_UBLAS_SINGULAR_CHECK
2047             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
2048 #else
2049             if (e1 () (n, n) == value_type/*zero*/())
2050                 singular ().raise ();
2051 #endif
2052             value_type t = e2 () (n) /= e1 () (n, n);
2053             if (t != value_type/*zero*/()) {
2054                 for (size_type m = n + 1; m < size; ++ m)
2055                     e2 () (m) -= e1 () (m, n) * t;
2056             }
2057         }
2058     }
2059     // Packed (proxy) case
2060     template<class E1, class E2>
2061     BOOST_UBLAS_INLINE
2062     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
2063                         lower_tag, column_major_tag, packed_proxy_tag) {
2064         typedef typename E2::size_type size_type;
2065         typedef typename E2::difference_type difference_type;
2066         typedef typename E2::value_type value_type;
2067
2068         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
2069         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
2070         size_type size = e2 ().size ();
2071         for (size_type n = 0; n < size; ++ n) {
2072 #ifndef BOOST_UBLAS_SINGULAR_CHECK
2073             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
2074 #else
2075             if (e1 () (n, n) == value_type/*zero*/())
2076                 singular ().raise ();
2077 #endif
2078             value_type t = e2 () (n) /= e1 () (n, n);
2079             if (t != value_type/*zero*/()) {
2080                 typename E1::const_iterator1 it1e1 (e1 ().find1 (1, n + 1, n));
2081                 typename E1::const_iterator1 it1e1_end (e1 ().find1 (1, e1 ().size1 (), n));
2082                 difference_type m (it1e1_end - it1e1);
2083                 while (-- m >= 0)
2084                     e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1;
2085             }
2086         }
2087     }
2088     // Sparse (proxy) case
2089     template<class E1, class E2>
2090     BOOST_UBLAS_INLINE
2091     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
2092                         lower_tag, column_major_tag, unknown_storage_tag) {
2093         typedef typename E2::size_type size_type;
2094         typedef typename E2::value_type value_type;
2095
2096         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
2097         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
2098         size_type size = e2 ().size ();
2099         for (size_type n = 0; n < size; ++ n) {
2100 #ifndef BOOST_UBLAS_SINGULAR_CHECK
2101             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
2102 #else
2103             if (e1 () (n, n) == value_type/*zero*/())
2104                 singular ().raise ();
2105 #endif
2106             value_type t = e2 () (n) /= e1 () (n, n);
2107             if (t != value_type/*zero*/()) {
2108                 typename E1::const_iterator1 it1e1 (e1 ().find1 (1, n + 1, n));
2109                 typename E1::const_iterator1 it1e1_end (e1 ().find1 (1, e1 ().size1 (), n));
2110                 while (it1e1 != it1e1_end)
2111                     e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1;
2112             }
2113         }
2114     }
2115
2116     // Dense (proxy) case
2117     template<class E1, class E2>
2118     BOOST_UBLAS_INLINE
2119     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
2120                         lower_tag, row_major_tag, dense_proxy_tag) {
2121         typedef typename E2::size_type size_type;
2122         typedef typename E2::value_type value_type;
2123
2124         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
2125         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
2126         size_type size = e2 ().size ();
2127         for (size_type n = 0; n < size; ++ n) {
2128 #ifndef BOOST_UBLAS_SINGULAR_CHECK
2129             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
2130 #else
2131             if (e1 () (n, n) == value_type/*zero*/())
2132                 singular ().raise ();
2133 #endif
2134             value_type t = e2 () (n) /= e1 () (n, n);
2135             if (t != value_type/*zero*/()) {
2136                 for (size_type m = n + 1; m < size; ++ m)
2137                     e2 () (m) -= e1 () (m, n) * t;
2138             }
2139         }
2140     }
2141     // Packed (proxy) case
2142     template<class E1, class E2>
2143     BOOST_UBLAS_INLINE
2144     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
2145                         lower_tag, row_major_tag, packed_proxy_tag) {
2146         typedef typename E2::size_type size_type;
2147         typedef typename E2::value_type value_type;
2148
2149         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
2150         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
2151         size_type size = e2 ().size ();
2152         for (size_type n = 0; n < size; ++ n) {
2153 #ifndef BOOST_UBLAS_SINGULAR_CHECK
2154             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
2155 #else
2156             if (e1 () (n, n) == value_type/*zero*/())
2157                 singular ().raise ();
2158 #endif
2159             value_type t = e2 () (n);
2160             typename E1::const_iterator2 it2e1 (e1 ().find2 (1, n, 0));
2161             typename E1::const_iterator2 it2e1_end (e1 ().find2 (1, n, n));
2162             while (it2e1 != it2e1_end) {
2163               t -= *it2e1 * e2 () (it2e1.index2());
2164               ++ it2e1;
2165             }
2166             e2() (n) = t / e1 () (n, n);
2167         }
2168     }
2169     // Sparse (proxy) case
2170     template<class E1, class E2>
2171     BOOST_UBLAS_INLINE
2172     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
2173                         lower_tag, row_major_tag, unknown_storage_tag) {
2174         typedef typename E2::size_type size_type;
2175         typedef typename E2::value_type value_type;
2176
2177         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
2178         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
2179         size_type size = e2 ().size ();
2180         for (size_type n = 0; n < size; ++ n) {
2181 #ifndef BOOST_UBLAS_SINGULAR_CHECK
2182             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
2183 #else
2184             if (e1 () (n, n) == value_type/*zero*/())
2185                 singular ().raise ();
2186 #endif
2187             value_type t = e2 () (n);
2188             typename E1::const_iterator2 it2e1 (e1 ().find2 (1, n, 0));
2189             typename E1::const_iterator2 it2e1_end (e1 ().find2 (1, n, n));
2190             while (it2e1 != it2e1_end) {
2191               t -= *it2e1 * e2 () (it2e1.index2());
2192               ++ it2e1;
2193             }
2194             e2() (n) = t / e1 () (n, n);
2195         }
2196     }
2197
2198     // Redirectors :-)
2199     template<class E1, class E2>
2200     BOOST_UBLAS_INLINE
2201     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
2202                         lower_tag, column_major_tag) {
2203         typedef typename E1::storage_category storage_category;
2204         inplace_solve (e1, e2,
2205                        lower_tag (), column_major_tag (), storage_category ());
2206     }
2207     template<class E1, class E2>
2208     BOOST_UBLAS_INLINE
2209     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
2210                         lower_tag, row_major_tag) {
2211         typedef typename E1::storage_category storage_category;
2212         inplace_solve (e1, e2,
2213                        lower_tag (), row_major_tag (), storage_category ());
2214     }
2215     // Dispatcher
2216     template<class E1, class E2>
2217     BOOST_UBLAS_INLINE
2218     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
2219                         lower_tag) {
2220         typedef typename E1::orientation_category orientation_category;
2221         inplace_solve (e1, e2,
2222                        lower_tag (), orientation_category ());
2223     }
2224     template<class E1, class E2>
2225     BOOST_UBLAS_INLINE
2226     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
2227                         unit_lower_tag) {
2228         typedef typename E1::orientation_category orientation_category;
2229         inplace_solve (triangular_adaptor<const E1, unit_lower> (e1 ()), e2,
2230                        unit_lower_tag (), orientation_category ());
2231     }
2232
2233     // Dense (proxy) case
2234     template<class E1, class E2>
2235     BOOST_UBLAS_INLINE
2236     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
2237                         upper_tag, column_major_tag, dense_proxy_tag) {
2238         typedef typename E2::size_type size_type;
2239         typedef typename E2::difference_type difference_type;
2240         typedef typename E2::value_type value_type;
2241
2242         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
2243         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
2244         size_type size = e2 ().size ();
2245         for (difference_type n = size - 1; n >= 0; -- n) {
2246 #ifndef BOOST_UBLAS_SINGULAR_CHECK
2247             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
2248 #else
2249             if (e1 () (n, n) == value_type/*zero*/())
2250                 singular ().raise ();
2251 #endif
2252             value_type t = e2 () (n) /= e1 () (n, n);
2253             if (t != value_type/*zero*/()) {
2254                 for (difference_type m = n - 1; m >= 0; -- m)
2255                     e2 () (m) -= e1 () (m, n) * t;
2256             }
2257         }
2258     }
2259     // Packed (proxy) case
2260     template<class E1, class E2>
2261     BOOST_UBLAS_INLINE
2262     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
2263                         upper_tag, column_major_tag, packed_proxy_tag) {
2264         typedef typename E2::size_type size_type;
2265         typedef typename E2::difference_type difference_type;
2266         typedef typename E2::value_type value_type;
2267
2268         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
2269         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
2270         size_type size = e2 ().size ();
2271         for (difference_type n = size - 1; n >= 0; -- n) {
2272 #ifndef BOOST_UBLAS_SINGULAR_CHECK
2273             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
2274 #else
2275             if (e1 () (n, n) == value_type/*zero*/())
2276                 singular ().raise ();
2277 #endif
2278             value_type t = e2 () (n) /= e1 () (n, n);
2279             if (t != value_type/*zero*/()) {
2280                 typename E1::const_reverse_iterator1 it1e1 (e1 ().find1 (1, n, n));
2281                 typename E1::const_reverse_iterator1 it1e1_rend (e1 ().find1 (1, 0, n));
2282                 while (it1e1 != it1e1_rend) {
2283                   e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1;
2284                 }
2285             }
2286         }
2287     }
2288     // Sparse (proxy) case
2289     template<class E1, class E2>
2290     BOOST_UBLAS_INLINE
2291     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
2292                         upper_tag, column_major_tag, unknown_storage_tag) {
2293         typedef typename E2::size_type size_type;
2294         typedef typename E2::difference_type difference_type;
2295         typedef typename E2::value_type value_type;
2296
2297         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
2298         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
2299         size_type size = e2 ().size ();
2300         for (difference_type n = size - 1; n >= 0; -- n) {
2301 #ifndef BOOST_UBLAS_SINGULAR_CHECK
2302             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
2303 #else
2304             if (e1 () (n, n) == value_type/*zero*/())
2305                 singular ().raise ();
2306 #endif
2307             value_type t = e2 () (n) /= e1 () (n, n);
2308             if (t != value_type/*zero*/()) {
2309                 typename E1::const_reverse_iterator1 it1e1 (e1 ().find1 (1, n, n));
2310                 typename E1::const_reverse_iterator1 it1e1_rend (e1 ().find1 (1, 0, n));
2311                 while (it1e1 != it1e1_rend) {
2312                   e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1;
2313                 }
2314             }
2315         }
2316     }
2317
2318     // Dense (proxy) case
2319     template<class E1, class E2>
2320     BOOST_UBLAS_INLINE
2321     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
2322                         upper_tag, row_major_tag, dense_proxy_tag) {
2323         typedef typename E2::size_type size_type;
2324         typedef typename E2::difference_type difference_type;
2325         typedef typename E2::value_type value_type;
2326
2327         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
2328         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
2329         size_type size = e1 ().size1 ();
2330         for (difference_type n = size-1; n >=0; -- n) {
2331 #ifndef BOOST_UBLAS_SINGULAR_CHECK
2332             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
2333 #else
2334             if (e1 () (n, n) == value_type/*zero*/())
2335                 singular ().raise ();
2336 #endif
2337             value_type t = e2 () (n);
2338             for (difference_type m = n + 1; m < static_cast<difference_type>(e1 ().size2()); ++ m) {
2339               t -= e1 () (n, m)  * e2 () (m);
2340             }
2341             e2() (n) = t / e1 () (n, n);
2342         }
2343     }
2344     // Packed (proxy) case
2345     template<class E1, class E2>
2346     BOOST_UBLAS_INLINE
2347     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
2348                         upper_tag, row_major_tag, packed_proxy_tag) {
2349         typedef typename E2::size_type size_type;
2350         typedef typename E2::difference_type difference_type;
2351         typedef typename E2::value_type value_type;
2352
2353         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
2354         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
2355         size_type size = e1 ().size1 ();
2356         for (difference_type n = size-1; n >=0; -- n) {
2357 #ifndef BOOST_UBLAS_SINGULAR_CHECK
2358             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
2359 #else
2360             if (e1 () (n, n) == value_type/*zero*/())
2361                 singular ().raise ();
2362 #endif
2363             value_type t = e2 () (n);
2364             typename E1::const_iterator2 it2e1 (e1 ().find2 (1, n, n+1));
2365             typename E1::const_iterator2 it2e1_end (e1 ().find2 (1, n, e1 ().size2 ()));
2366             while (it2e1 != it2e1_end) {
2367               t -= *it2e1 * e2 () (it2e1.index2());
2368               ++ it2e1;
2369             }
2370             e2() (n) = t / e1 () (n, n);
2371
2372         }
2373     }
2374     // Sparse (proxy) case
2375     template<class E1, class E2>
2376     BOOST_UBLAS_INLINE
2377     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
2378                         upper_tag, row_major_tag, unknown_storage_tag) {
2379         typedef typename E2::size_type size_type;
2380         typedef typename E2::difference_type difference_type;
2381         typedef typename E2::value_type value_type;
2382
2383         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
2384         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
2385         size_type size = e1 ().size1 ();
2386         for (difference_type n = size-1; n >=0; -- n) {
2387 #ifndef BOOST_UBLAS_SINGULAR_CHECK
2388             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
2389 #else
2390             if (e1 () (n, n) == value_type/*zero*/())
2391                 singular ().raise ();
2392 #endif
2393             value_type t = e2 () (n);
2394             typename E1::const_iterator2 it2e1 (e1 ().find2 (1, n, n+1));
2395             typename E1::const_iterator2 it2e1_end (e1 ().find2 (1, n, e1 ().size2 ()));
2396             while (it2e1 != it2e1_end) {
2397               t -= *it2e1 * e2 () (it2e1.index2());
2398               ++ it2e1;
2399             }
2400             e2() (n) = t / e1 () (n, n);
2401
2402         }
2403     }
2404
2405     // Redirectors :-)
2406     template<class E1, class E2>
2407     BOOST_UBLAS_INLINE
2408     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
2409                         upper_tag, column_major_tag) {
2410         typedef typename E1::storage_category storage_category;
2411         inplace_solve (e1, e2,
2412                        upper_tag (), column_major_tag (), storage_category ());
2413     }
2414     template<class E1, class E2>
2415     BOOST_UBLAS_INLINE
2416     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
2417                         upper_tag, row_major_tag) {
2418         typedef typename E1::storage_category storage_category;
2419         inplace_solve (e1, e2,
2420                        upper_tag (), row_major_tag (), storage_category ());
2421     }
2422     // Dispatcher
2423     template<class E1, class E2>
2424     BOOST_UBLAS_INLINE
2425     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
2426                         upper_tag) {
2427         typedef typename E1::orientation_category orientation_category;
2428         inplace_solve (e1, e2,
2429                        upper_tag (), orientation_category ());
2430     }
2431     template<class E1, class E2>
2432     BOOST_UBLAS_INLINE
2433     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
2434                         unit_upper_tag) {
2435         typedef typename E1::orientation_category orientation_category;
2436         inplace_solve (triangular_adaptor<const E1, unit_upper> (e1 ()), e2,
2437                        unit_upper_tag (), orientation_category ());
2438     }
2439
2440     template<class E1, class E2, class C>
2441     BOOST_UBLAS_INLINE
2442     typename matrix_vector_solve_traits<E1, E2>::result_type
2443     solve (const matrix_expression<E1> &e1,
2444            const vector_expression<E2> &e2,
2445            C) {
2446         typename matrix_vector_solve_traits<E1, E2>::result_type r (e2);
2447         inplace_solve (e1, r, C ());
2448         return r;
2449     }
2450
2451
2452     // Redirectors :-)
2453     template<class E1, class E2>
2454     BOOST_UBLAS_INLINE
2455     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
2456                         lower_tag, row_major_tag) {
2457         typedef typename E2::storage_category storage_category;
2458         inplace_solve (trans(e2), e1,
2459                        upper_tag (), column_major_tag (), storage_category ());
2460     }
2461     template<class E1, class E2>
2462     BOOST_UBLAS_INLINE
2463     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
2464                         lower_tag, column_major_tag) {
2465         typedef typename E2::storage_category storage_category;
2466         inplace_solve (trans (e2), e1,
2467                        upper_tag (), row_major_tag (), storage_category ());
2468     }
2469     // Dispatcher
2470     template<class E1, class E2>
2471     BOOST_UBLAS_INLINE
2472     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
2473                         lower_tag) {
2474         typedef typename E2::orientation_category orientation_category;
2475         inplace_solve (e1, e2,
2476                        lower_tag (), orientation_category ());
2477     }
2478     template<class E1, class E2>
2479     BOOST_UBLAS_INLINE
2480     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
2481                         unit_lower_tag) {
2482         typedef typename E2::orientation_category orientation_category;
2483         inplace_solve (e1, triangular_adaptor<const E2, unit_lower> (e2 ()),
2484                        unit_lower_tag (), orientation_category ());
2485     }
2486
2487
2488     // Redirectors :-)
2489     template<class E1, class E2>
2490     BOOST_UBLAS_INLINE
2491     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
2492                         upper_tag, row_major_tag) {
2493         typedef typename E2::storage_category storage_category;
2494         inplace_solve (trans(e2), e1,
2495                        lower_tag (), column_major_tag (), storage_category ());
2496     }
2497     template<class E1, class E2>
2498     BOOST_UBLAS_INLINE
2499     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
2500                         upper_tag, column_major_tag) {
2501         typedef typename E2::storage_category storage_category;
2502         inplace_solve (trans (e2), e1,
2503                        lower_tag (), row_major_tag (), storage_category ());
2504     }
2505     // Dispatcher
2506     template<class E1, class E2>
2507     BOOST_UBLAS_INLINE
2508     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
2509                         upper_tag) {
2510         typedef typename E2::orientation_category orientation_category;
2511         inplace_solve (e1, e2,
2512                        upper_tag (), orientation_category ());
2513     }
2514     template<class E1, class E2>
2515     BOOST_UBLAS_INLINE
2516     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
2517                         unit_upper_tag) {
2518         typedef typename E2::orientation_category orientation_category;
2519         inplace_solve (e1, triangular_adaptor<const E2, unit_upper> (e2 ()),
2520                        unit_upper_tag (), orientation_category ());
2521     }
2522
2523     template<class E1, class E2, class C>
2524     BOOST_UBLAS_INLINE
2525     typename matrix_vector_solve_traits<E1, E2>::result_type
2526     solve (const vector_expression<E1> &e1,
2527            const matrix_expression<E2> &e2,
2528            C) {
2529         typename matrix_vector_solve_traits<E1, E2>::result_type r (e1);
2530         inplace_solve (r, e2, C ());
2531         return r;
2532     }
2533
2534     template<class E1, class E2>
2535     struct matrix_matrix_solve_traits {
2536         typedef typename promote_traits<typename E1::value_type, typename E2::value_type>::promote_type promote_type;
2537         typedef matrix<promote_type> result_type;
2538     };
2539
2540     // Operations:
2541     //  k * n * (n - 1) / 2 + k * n = k * n * (n + 1) / 2 multiplications,
2542     //  k * n * (n - 1) / 2 additions
2543
2544     // Dense (proxy) case
2545     template<class E1, class E2>
2546     BOOST_UBLAS_INLINE
2547     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
2548                         lower_tag, dense_proxy_tag) {
2549         typedef typename E2::size_type size_type;
2550         typedef typename E2::value_type value_type;
2551
2552         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
2553         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());
2554         size_type size1 = e2 ().size1 ();
2555         size_type size2 = e2 ().size2 ();
2556         for (size_type n = 0; n < size1; ++ n) {
2557 #ifndef BOOST_UBLAS_SINGULAR_CHECK
2558             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
2559 #else
2560             if (e1 () (n, n) == value_type/*zero*/())
2561                 singular ().raise ();
2562 #endif
2563             for (size_type l = 0; l < size2; ++ l) {
2564                 value_type t = e2 () (n, l) /= e1 () (n, n);
2565                 if (t != value_type/*zero*/()) {
2566                     for (size_type m = n + 1; m < size1; ++ m)
2567                         e2 () (m, l) -= e1 () (m, n) * t;
2568                 }
2569             }
2570         }
2571     }
2572     // Packed (proxy) case
2573     template<class E1, class E2>
2574     BOOST_UBLAS_INLINE
2575     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
2576                         lower_tag, packed_proxy_tag) {
2577         typedef typename E2::size_type size_type;
2578         typedef typename E2::difference_type difference_type;
2579         typedef typename E2::value_type value_type;
2580
2581         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
2582         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());
2583         size_type size1 = e2 ().size1 ();
2584         size_type size2 = e2 ().size2 ();
2585         for (size_type n = 0; n < size1; ++ n) {
2586 #ifndef BOOST_UBLAS_SINGULAR_CHECK
2587             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
2588 #else
2589             if (e1 () (n, n) == value_type/*zero*/())
2590                 singular ().raise ();
2591 #endif
2592             for (size_type l = 0; l < size2; ++ l) {
2593                 value_type t = e2 () (n, l) /= e1 () (n, n);
2594                 if (t != value_type/*zero*/()) {
2595                     typename E1::const_iterator1 it1e1 (e1 ().find1 (1, n + 1, n));
2596                     typename E1::const_iterator1 it1e1_end (e1 ().find1 (1, e1 ().size1 (), n));
2597                     difference_type m (it1e1_end - it1e1);
2598                     while (-- m >= 0)
2599                         e2 () (it1e1.index1 (), l) -= *it1e1 * t, ++ it1e1;
2600                 }
2601             }
2602         }
2603     }
2604     // Sparse (proxy) case
2605     template<class E1, class E2>
2606     BOOST_UBLAS_INLINE
2607     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
2608                         lower_tag, unknown_storage_tag) {
2609         typedef typename E2::size_type size_type;
2610         typedef typename E2::value_type value_type;
2611
2612         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
2613         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());
2614         size_type size1 = e2 ().size1 ();
2615         size_type size2 = e2 ().size2 ();
2616         for (size_type n = 0; n < size1; ++ n) {
2617 #ifndef BOOST_UBLAS_SINGULAR_CHECK
2618             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
2619 #else
2620             if (e1 () (n, n) == value_type/*zero*/())
2621                 singular ().raise ();
2622 #endif
2623             for (size_type l = 0; l < size2; ++ l) {
2624                 value_type t = e2 () (n, l) /= e1 () (n, n);
2625                 if (t != value_type/*zero*/()) {
2626                     typename E1::const_iterator1 it1e1 (e1 ().find1 (1, n + 1, n));
2627                     typename E1::const_iterator1 it1e1_end (e1 ().find1 (1, e1 ().size1 (), n));
2628                     while (it1e1 != it1e1_end)
2629                         e2 () (it1e1.index1 (), l) -= *it1e1 * t, ++ it1e1;
2630                 }
2631             }
2632         }
2633     }
2634     // Dispatcher
2635     template<class E1, class E2>
2636     BOOST_UBLAS_INLINE
2637     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
2638                         lower_tag) {
2639         typedef typename E1::storage_category dispatch_category;
2640         inplace_solve (e1, e2,
2641                        lower_tag (), dispatch_category ());
2642     }
2643     template<class E1, class E2>
2644     BOOST_UBLAS_INLINE
2645     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
2646                         unit_lower_tag) {
2647         typedef typename E1::storage_category dispatch_category;
2648         inplace_solve (triangular_adaptor<const E1, unit_lower> (e1 ()), e2,
2649                        unit_lower_tag (), dispatch_category ());
2650     }
2651
2652     // Dense (proxy) case
2653     template<class E1, class E2>
2654     BOOST_UBLAS_INLINE
2655     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
2656                         upper_tag, dense_proxy_tag) {
2657         typedef typename E2::size_type size_type;
2658         typedef typename E2::difference_type difference_type;
2659         typedef typename E2::value_type value_type;
2660
2661         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
2662         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());
2663         size_type size1 = e2 ().size1 ();
2664         size_type size2 = e2 ().size2 ();
2665         for (difference_type n = size1 - 1; n >= 0; -- n) {
2666 #ifndef BOOST_UBLAS_SINGULAR_CHECK
2667             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
2668 #else
2669             if (e1 () (n, n) == value_type/*zero*/())
2670                 singular ().raise ();
2671 #endif
2672             for (difference_type l = size2 - 1; l >= 0; -- l) {
2673                 value_type t = e2 () (n, l) /= e1 () (n, n);
2674                 if (t != value_type/*zero*/()) {
2675                     for (difference_type m = n - 1; m >= 0; -- m)
2676                         e2 () (m, l) -= e1 () (m, n) * t;
2677                 }
2678             }
2679         }
2680     }
2681     // Packed (proxy) case
2682     template<class E1, class E2>
2683     BOOST_UBLAS_INLINE
2684     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
2685                         upper_tag, packed_proxy_tag) {
2686         typedef typename E2::size_type size_type;
2687         typedef typename E2::difference_type difference_type;
2688         typedef typename E2::value_type value_type;
2689
2690         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
2691         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());
2692         size_type size1 = e2 ().size1 ();
2693         size_type size2 = e2 ().size2 ();
2694         for (difference_type n = size1 - 1; n >= 0; -- n) {
2695 #ifndef BOOST_UBLAS_SINGULAR_CHECK
2696             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
2697 #else
2698             if (e1 () (n, n) == value_type/*zero*/())
2699                 singular ().raise ();
2700 #endif
2701             for (difference_type l = size2 - 1; l >= 0; -- l) {
2702                 value_type t = e2 () (n, l) /= e1 () (n, n);
2703                 if (t != value_type/*zero*/()) {
2704                     typename E1::const_reverse_iterator1 it1e1 (e1 ().find1 (1, n, n));
2705                     typename E1::const_reverse_iterator1 it1e1_rend (e1 ().find1 (1, 0, n));
2706                     difference_type m (it1e1_rend - it1e1);
2707                     while (-- m >= 0)
2708                         e2 () (it1e1.index1 (), l) -= *it1e1 * t, ++ it1e1;
2709                 }
2710             }
2711         }
2712     }
2713     // Sparse (proxy) case
2714     template<class E1, class E2>
2715     BOOST_UBLAS_INLINE
2716     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
2717                         upper_tag, unknown_storage_tag) {
2718         typedef typename E2::size_type size_type;
2719         typedef typename E2::difference_type difference_type;
2720         typedef typename E2::value_type value_type;
2721
2722         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
2723         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());
2724         size_type size1 = e2 ().size1 ();
2725         size_type size2 = e2 ().size2 ();
2726         for (difference_type n = size1 - 1; n >= 0; -- n) {
2727 #ifndef BOOST_UBLAS_SINGULAR_CHECK
2728             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
2729 #else
2730             if (e1 () (n, n) == value_type/*zero*/())
2731                 singular ().raise ();
2732 #endif
2733             for (difference_type l = size2 - 1; l >= 0; -- l) {
2734                 value_type t = e2 () (n, l) /= e1 () (n, n);
2735                 if (t != value_type/*zero*/()) {
2736                     typename E1::const_reverse_iterator1 it1e1 (e1 ().find1 (1, n, n));
2737                     typename E1::const_reverse_iterator1 it1e1_rend (e1 ().find1 (1, 0, n));
2738                     while (it1e1 != it1e1_rend)
2739                         e2 () (it1e1.index1 (), l) -= *it1e1 * t, ++ it1e1;
2740                 }
2741             }
2742         }
2743     }
2744     // Dispatcher
2745     template<class E1, class E2>
2746     BOOST_UBLAS_INLINE
2747     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
2748                         upper_tag) {
2749         typedef typename E1::storage_category dispatch_category;
2750         inplace_solve (e1, e2,
2751                        upper_tag (), dispatch_category ());
2752     }
2753     template<class E1, class E2>
2754     BOOST_UBLAS_INLINE
2755     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
2756                         unit_upper_tag) {
2757         typedef typename E1::storage_category dispatch_category;
2758         inplace_solve (triangular_adaptor<const E1, unit_upper> (e1 ()), e2,
2759                        unit_upper_tag (), dispatch_category ());
2760     }
2761
2762     template<class E1, class E2, class C>
2763     BOOST_UBLAS_INLINE
2764     typename matrix_matrix_solve_traits<E1, E2>::result_type
2765     solve (const matrix_expression<E1> &e1,
2766            const matrix_expression<E2> &e2,
2767            C) {
2768         typename matrix_matrix_solve_traits<E1, E2>::result_type r (e2);
2769         inplace_solve (e1, r, C ());
2770         return r;
2771     }
2772
2773 }}}
2774
2775 #endif