Symaptic: os/ossrv/ossrv_pub/boost_apis/boost/numeric/ublas/triangular.hpp@260cb5ec6c19

     1 //

     2 //  Copyright (c) 2000-2002

     3 //  Joerg Walter, Mathias Koch

     4 //

     5 //  Permission to use, copy, modify, distribute and sell this software

     6 //  and its documentation for any purpose is hereby granted without fee,

     7 //  provided that the above copyright notice appear in all copies and

     8 //  that both that copyright notice and this permission notice appear

     9 //  in supporting documentation.  The authors make no representations

    10 //  about the suitability of this software for any purpose.

    11 //  It is provided "as is" without express or implied warranty.

    12 //

    13 //  The authors gratefully acknowledge the support of

    14 //  GeNeSys mbH & Co. KG in producing this work.

    15 //

    17 #ifndef _BOOST_UBLAS_TRIANGULAR_

    18 #define _BOOST_UBLAS_TRIANGULAR_

    20 #include <boost/numeric/ublas/matrix.hpp>

    21 #include <boost/numeric/ublas/detail/temporary.hpp>

    22 #include <boost/type_traits/remove_const.hpp>

    24 // Iterators based on ideas of Jeremy Siek

    26 namespace boost { namespace numeric { namespace ublas {

    28     // Array based triangular matrix class

    29     template<class T, class TRI, class L, class A>

    30     class triangular_matrix:

    31         public matrix_container<triangular_matrix<T, TRI, L, A> > {

    33         typedef T *pointer;

    34         typedef TRI triangular_type;

    35         typedef L layout_type;

    36         typedef triangular_matrix<T, TRI, L, A> self_type;

    37     public:

    38 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS

    39         using matrix_container<self_type>::operator ();

    40 #endif

    41         typedef typename A::size_type size_type;

    42         typedef typename A::difference_type difference_type;

    43         typedef T value_type;

    44         typedef const T &const_reference;

    45         typedef T &reference;

    46         typedef A array_type;

    48         typedef const matrix_reference<const self_type> const_closure_type;

    49         typedef matrix_reference<self_type> closure_type;

    50         typedef vector<T, A> vector_temporary_type;

    51         typedef matrix<T, L, A> matrix_temporary_type;  // general sub-matrix

    52         typedef packed_tag storage_category;

    53         typedef typename L::orientation_category orientation_category;

    55         // Construction and destruction

    56         BOOST_UBLAS_INLINE

    57         triangular_matrix ():

    58             matrix_container<self_type> (),

    59             size1_ (0), size2_ (0), data_ (0) {}

    60         BOOST_UBLAS_INLINE

    61         triangular_matrix (size_type size1, size_type size2):

    62             matrix_container<self_type> (),

    63             size1_ (size1), size2_ (size2), data_ (triangular_type::packed_size (layout_type (), size1, size2)) {

    64         }

    65         BOOST_UBLAS_INLINE

    66         triangular_matrix (size_type size1, size_type size2, const array_type &data):

    67             matrix_container<self_type> (),

    68             size1_ (size1), size2_ (size2), data_ (data) {}

    69         BOOST_UBLAS_INLINE

    70         triangular_matrix (const triangular_matrix &m):

    71             matrix_container<self_type> (),

    72             size1_ (m.size1_), size2_ (m.size2_), data_ (m.data_) {}

    73         template<class AE>

    74         BOOST_UBLAS_INLINE

    75         triangular_matrix (const matrix_expression<AE> &ae):

    76             matrix_container<self_type> (),

    77             size1_ (ae ().size1 ()), size2_ (ae ().size2 ()),

    78             data_ (triangular_type::packed_size (layout_type (), size1_, size2_)) {

    79             matrix_assign<scalar_assign> (*this, ae);

    80         }

    82         // Accessors

    83         BOOST_UBLAS_INLINE

    84         size_type size1 () const {

    85             return size1_;

    86         }

    87         BOOST_UBLAS_INLINE

    88         size_type size2 () const {

    89             return size2_;

    90         }

    92         // Storage accessors

    93         BOOST_UBLAS_INLINE

    94         const array_type &data () const {

    95             return data_;

    96         }

    97         BOOST_UBLAS_INLINE

    98         array_type &data () {

    99             return data_;

   100         }

   102         // Resizing

   103         BOOST_UBLAS_INLINE

   104         void resize (size_type size1, size_type size2, bool preserve = true) {

   105             if (preserve) {

   106                 self_type temporary (size1, size2);

   107                 detail::matrix_resize_preserve<layout_type> (*this, temporary);

   108             }

   109             else {

   110                 data ().resize (triangular_type::packed_size (layout_type (), size1, size2));

   111                 size1_ = size1;

   112                 size2_ = size2;

   113             }

   114         }

   115         BOOST_UBLAS_INLINE

   116         void resize_packed_preserve (size_type size1, size_type size2) {

   117             size1_ = size1;

   118             size2_ = size2;

   119             data ().resize (triangular_type::packed_size (layout_type (), size1_, size2_), value_type ());

   120         }

   122         // Element access

   123         BOOST_UBLAS_INLINE

   124         const_reference operator () (size_type i, size_type j) const {

   125             BOOST_UBLAS_CHECK (i < size1_, bad_index ());

   126             BOOST_UBLAS_CHECK (j < size2_, bad_index ());

   127             if (triangular_type::other (i, j))

   128                 return data () [triangular_type::element (layout_type (), i, size1_, j, size2_)];

   129             else if (triangular_type::one (i, j))

   130                 return one_;

   131             else

   132                 return zero_;

   133         }

   134         BOOST_UBLAS_INLINE

   135         reference at_element (size_type i, size_type j) {

   136             BOOST_UBLAS_CHECK (i < size1_, bad_index ());

   137             BOOST_UBLAS_CHECK (j < size2_, bad_index ());

   138             return data () [triangular_type::element (layout_type (), i, size1_, j, size2_)];

   139         }

   140         BOOST_UBLAS_INLINE

   141         reference operator () (size_type i, size_type j) {

   142             BOOST_UBLAS_CHECK (i < size1_, bad_index ());

   143             BOOST_UBLAS_CHECK (j < size2_, bad_index ());

   144             if (triangular_type::other (i, j))

   145                 return data () [triangular_type::element (layout_type (), i, size1_, j, size2_)];

   146             else {

   147                 bad_index ().raise ();

   148                 // arbitary return value

   149                 return const_cast<reference>(zero_);

   150             }

   151         }

   153         // Element assignment

   154         BOOST_UBLAS_INLINE

   155         reference insert_element (size_type i, size_type j, const_reference t) {

   156             return (operator () (i, j) = t);

   157         }

   158         BOOST_UBLAS_INLINE

   159         void erase_element (size_type i, size_type j) {

   160             operator () (i, j) = value_type/*zero*/();

   161         }

   163         // Zeroing

   164         BOOST_UBLAS_INLINE

   165         void clear () {

   166             // data ().clear ();

   167             std::fill (data ().begin (), data ().end (), value_type/*zero*/());

   168         }

   170         // Assignment

   171         BOOST_UBLAS_INLINE

   172         triangular_matrix &operator = (const triangular_matrix &m) {

   173             size1_ = m.size1_;

   174             size2_ = m.size2_;

   175             data () = m.data ();

   176             return *this;

   177         }

   178         BOOST_UBLAS_INLINE

   179         triangular_matrix &assign_temporary (triangular_matrix &m) {

   180             swap (m);

   181             return *this;

   182         }

   183         template<class AE>

   184         BOOST_UBLAS_INLINE

   185         triangular_matrix &operator = (const matrix_expression<AE> &ae) {

   186             self_type temporary (ae);

   187             return assign_temporary (temporary);

   188         }

   189         template<class AE>

   190         BOOST_UBLAS_INLINE

   191         triangular_matrix &assign (const matrix_expression<AE> &ae) {

   192             matrix_assign<scalar_assign> (*this, ae);

   193             return *this;

   194         }

   195         template<class AE>

   196         BOOST_UBLAS_INLINE

   197         triangular_matrix& operator += (const matrix_expression<AE> &ae) {

   198             self_type temporary (*this + ae);

   199             return assign_temporary (temporary);

   200         }

   201         template<class AE>

   202         BOOST_UBLAS_INLINE

   203         triangular_matrix &plus_assign (const matrix_expression<AE> &ae) {

   204             matrix_assign<scalar_plus_assign> (*this, ae);

   205             return *this;

   206         }

   207         template<class AE>

   208         BOOST_UBLAS_INLINE

   209         triangular_matrix& operator -= (const matrix_expression<AE> &ae) {

   210             self_type temporary (*this - ae);

   211             return assign_temporary (temporary);

   212         }

   213         template<class AE>

   214         BOOST_UBLAS_INLINE

   215         triangular_matrix &minus_assign (const matrix_expression<AE> &ae) {

   216             matrix_assign<scalar_minus_assign> (*this, ae);

   217             return *this;

   218         }

   219         template<class AT>

   220         BOOST_UBLAS_INLINE

   221         triangular_matrix& operator *= (const AT &at) {

   222             matrix_assign_scalar<scalar_multiplies_assign> (*this, at);

   223             return *this;

   224         }

   225         template<class AT>

   226         BOOST_UBLAS_INLINE

   227         triangular_matrix& operator /= (const AT &at) {

   228             matrix_assign_scalar<scalar_divides_assign> (*this, at);

   229             return *this;

   230         }

   232         // Swapping

   233         BOOST_UBLAS_INLINE

   234         void swap (triangular_matrix &m) {

   235             if (this != &m) {

   236                 // BOOST_UBLAS_CHECK (size2_ == m.size2_, bad_size ());

   237                 std::swap (size1_, m.size1_);

   238                 std::swap (size2_, m.size2_);

   239                 data ().swap (m.data ());

   240             }

   241         }

   242         BOOST_UBLAS_INLINE

   243         friend void swap (triangular_matrix &m1, triangular_matrix &m2) {

   244             m1.swap (m2);

   245         }

   247         // Iterator types

   248 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR

   249         typedef indexed_iterator1<self_type, packed_random_access_iterator_tag> iterator1;

   250         typedef indexed_iterator2<self_type, packed_random_access_iterator_tag> iterator2;

   251         typedef indexed_const_iterator1<self_type, packed_random_access_iterator_tag> const_iterator1;

   252         typedef indexed_const_iterator2<self_type, packed_random_access_iterator_tag> const_iterator2;

   253 #else

   254         class const_iterator1;

   255         class iterator1;

   256         class const_iterator2;

   257         class iterator2;

   258 #endif

   259         typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;

   260         typedef reverse_iterator_base1<iterator1> reverse_iterator1;

   261         typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;

   262         typedef reverse_iterator_base2<iterator2> reverse_iterator2;

   264         // Element lookup

   265         BOOST_UBLAS_INLINE

   266         const_iterator1 find1 (int rank, size_type i, size_type j) const {

   267             if (rank == 1)

   268                 i = triangular_type::restrict1 (i, j);

   269             return const_iterator1 (*this, i, j);

   270         }

   271         BOOST_UBLAS_INLINE

   272         iterator1 find1 (int rank, size_type i, size_type j) {

   273             if (rank == 1)

   274                 i = triangular_type::mutable_restrict1 (i, j);

   275             return iterator1 (*this, i, j);

   276         }

   277         BOOST_UBLAS_INLINE

   278         const_iterator2 find2 (int rank, size_type i, size_type j) const {

   279             if (rank == 1)

   280                 j = triangular_type::restrict2 (i, j);

   281             return const_iterator2 (*this, i, j);

   282         }

   283         BOOST_UBLAS_INLINE

   284         iterator2 find2 (int rank, size_type i, size_type j) {

   285             if (rank == 1)

   286                 j = triangular_type::mutable_restrict2 (i, j);

   287             return iterator2 (*this, i, j);

   288         }

   290         // Iterators simply are indices.

   292 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR

   293         class const_iterator1:

   294             public container_const_reference<triangular_matrix>,

   295             public random_access_iterator_base<packed_random_access_iterator_tag,

   296                                                const_iterator1, value_type> {

   297         public:

   298             typedef typename triangular_matrix::value_type value_type;

   299             typedef typename triangular_matrix::difference_type difference_type;

   300             typedef typename triangular_matrix::const_reference reference;

   301             typedef const typename triangular_matrix::pointer pointer;

   303             typedef const_iterator2 dual_iterator_type;

   304             typedef const_reverse_iterator2 dual_reverse_iterator_type;

   306             // Construction and destruction

   307             BOOST_UBLAS_INLINE

   308             const_iterator1 ():

   309                 container_const_reference<self_type> (), it1_ (), it2_ () {}

   310             BOOST_UBLAS_INLINE

   311             const_iterator1 (const self_type &m, size_type it1, size_type it2):

   312                 container_const_reference<self_type> (m), it1_ (it1), it2_ (it2) {}

   313             BOOST_UBLAS_INLINE

   314             const_iterator1 (const iterator1 &it):

   315                 container_const_reference<self_type> (it ()), it1_ (it.it1_), it2_ (it.it2_) {}

   317             // Arithmetic

   318             BOOST_UBLAS_INLINE

   319             const_iterator1 &operator ++ () {

   320                 ++ it1_;

   321                 return *this;

   322             }

   323             BOOST_UBLAS_INLINE

   324             const_iterator1 &operator -- () {

   325                 -- it1_;

   326                 return *this;

   327             }

   328             BOOST_UBLAS_INLINE

   329             const_iterator1 &operator += (difference_type n) {

   330                 it1_ += n;

   331                 return *this;

   332             }

   333             BOOST_UBLAS_INLINE

   334             const_iterator1 &operator -= (difference_type n) {

   335                 it1_ -= n;

   336                 return *this;

   337             }

   338             BOOST_UBLAS_INLINE

   339             difference_type operator - (const const_iterator1 &it) const {

   340                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

   341                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());

   342                 return it1_ - it.it1_;

   343             }

   345             // Dereference

   346             BOOST_UBLAS_INLINE

   347             const_reference operator * () const {

   348                 return (*this) () (it1_, it2_);

   349             }

   350             BOOST_UBLAS_INLINE

   351             const_reference operator [] (difference_type n) const {

   352                 return *(*this + n);

   353             }

   355 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION

   356             BOOST_UBLAS_INLINE

   357 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

   358             typename self_type::

   359 #endif

   360             const_iterator2 begin () const {

   361                 return (*this) ().find2 (1, it1_, 0);

   362             }

   363             BOOST_UBLAS_INLINE

   364 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

   365             typename self_type::

   366 #endif

   367             const_iterator2 end () const {

   368                 return (*this) ().find2 (1, it1_, (*this) ().size2 ());

   369             }

   370             BOOST_UBLAS_INLINE

   371 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

   372             typename self_type::

   373 #endif

   374             const_reverse_iterator2 rbegin () const {

   375                 return const_reverse_iterator2 (end ());

   376             }

   377             BOOST_UBLAS_INLINE

   378 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

   379             typename self_type::

   380 #endif

   381             const_reverse_iterator2 rend () const {

   382                 return const_reverse_iterator2 (begin ());

   383             }

   384 #endif

   386             // Indices

   387             BOOST_UBLAS_INLINE

   388             size_type index1 () const {

   389                 return it1_;

   390             }

   391             BOOST_UBLAS_INLINE

   392             size_type index2 () const {

   393                 return it2_;

   394             }

   396             // Assignment

   397             BOOST_UBLAS_INLINE

   398             const_iterator1 &operator = (const const_iterator1 &it) {

   399                 container_const_reference<self_type>::assign (&it ());

   400                 it1_ = it.it1_;

   401                 it2_ = it.it2_;

   402                 return *this;

   403             }

   405             // Comparison

   406             BOOST_UBLAS_INLINE

   407             bool operator == (const const_iterator1 &it) const {

   408                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

   409                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());

   410                 return it1_ == it.it1_;

   411             }

   412             BOOST_UBLAS_INLINE

   413             bool operator < (const const_iterator1 &it) const {

   414                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

   415                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());

   416                 return it1_ < it.it1_;

   417             }

   419         private:

   420             size_type it1_;

   421             size_type it2_;

   422         };

   423 #endif

   425         BOOST_UBLAS_INLINE

   426         const_iterator1 begin1 () const {

   427             return find1 (0, 0, 0);

   428         }

   429         BOOST_UBLAS_INLINE

   430         const_iterator1 end1 () const {

   431             return find1 (0, size1_, 0);

   432         }

   434 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR

   435         class iterator1:

   436             public container_reference<triangular_matrix>,

   437             public random_access_iterator_base<packed_random_access_iterator_tag,

   438                                                iterator1, value_type> {

   439         public:

   440             typedef typename triangular_matrix::value_type value_type;

   441             typedef typename triangular_matrix::difference_type difference_type;

   442             typedef typename triangular_matrix::reference reference;

   443             typedef typename triangular_matrix::pointer pointer;

   445             typedef iterator2 dual_iterator_type;

   446             typedef reverse_iterator2 dual_reverse_iterator_type;

   448             // Construction and destruction

   449             BOOST_UBLAS_INLINE

   450             iterator1 ():

   451                 container_reference<self_type> (), it1_ (), it2_ () {}

   452             BOOST_UBLAS_INLINE

   453             iterator1 (self_type &m, size_type it1, size_type it2):

   454                 container_reference<self_type> (m), it1_ (it1), it2_ (it2) {}

   456             // Arithmetic

   457             BOOST_UBLAS_INLINE

   458             iterator1 &operator ++ () {

   459                 ++ it1_;

   460                 return *this;

   461             }

   462             BOOST_UBLAS_INLINE

   463             iterator1 &operator -- () {

   464                 -- it1_;

   465                 return *this;

   466             }

   467             BOOST_UBLAS_INLINE

   468             iterator1 &operator += (difference_type n) {

   469                 it1_ += n;

   470                 return *this;

   471             }

   472             BOOST_UBLAS_INLINE

   473             iterator1 &operator -= (difference_type n) {

   474                 it1_ -= n;

   475                 return *this;

   476             }

   477             BOOST_UBLAS_INLINE

   478             difference_type operator - (const iterator1 &it) const {

   479                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

   480                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());

   481                 return it1_ - it.it1_;

   482             }

   484             // Dereference

   485             BOOST_UBLAS_INLINE

   486             reference operator * () const {

   487                 return (*this) () (it1_, it2_);

   488             }

   489             BOOST_UBLAS_INLINE

   490             reference operator [] (difference_type n) const {

   491                 return *(*this + n);

   492             }

   494 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION

   495             BOOST_UBLAS_INLINE

   496 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

   497             typename self_type::

   498 #endif

   499             iterator2 begin () const {

   500                 return (*this) ().find2 (1, it1_, 0);

   501             }

   502             BOOST_UBLAS_INLINE

   503 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

   504             typename self_type::

   505 #endif

   506             iterator2 end () const {

   507                 return (*this) ().find2 (1, it1_, (*this) ().size2 ());

   508             }

   509             BOOST_UBLAS_INLINE

   510 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

   511             typename self_type::

   512 #endif

   513             reverse_iterator2 rbegin () const {

   514                 return reverse_iterator2 (end ());

   515             }

   516             BOOST_UBLAS_INLINE

   517 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

   518             typename self_type::

   519 #endif

   520             reverse_iterator2 rend () const {

   521                 return reverse_iterator2 (begin ());

   522             }

   523 #endif

   525             // Indices

   526             BOOST_UBLAS_INLINE

   527             size_type index1 () const {

   528                 return it1_;

   529             }

   530             BOOST_UBLAS_INLINE

   531             size_type index2 () const {

   532                 return it2_;

   533             }

   535             // Assignment

   536             BOOST_UBLAS_INLINE

   537             iterator1 &operator = (const iterator1 &it) {

   538                 container_reference<self_type>::assign (&it ());

   539                 it1_ = it.it1_;

   540                 it2_ = it.it2_;

   541                 return *this;

   542             }

   544             // Comparison

   545             BOOST_UBLAS_INLINE

   546             bool operator == (const iterator1 &it) const {

   547                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

   548                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());

   549                 return it1_ == it.it1_;

   550             }

   551             BOOST_UBLAS_INLINE

   552             bool operator < (const iterator1 &it) const {

   553                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

   554                 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());

   555                 return it1_ < it.it1_;

   556             }

   558         private:

   559             size_type it1_;

   560             size_type it2_;

   562             friend class const_iterator1;

   563         };

   564 #endif

   566         BOOST_UBLAS_INLINE

   567         iterator1 begin1 () {

   568             return find1 (0, 0, 0);

   569         }

   570         BOOST_UBLAS_INLINE

   571         iterator1 end1 () {

   572             return find1 (0, size1_, 0);

   573         }

   575 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR

   576         class const_iterator2:

   577             public container_const_reference<triangular_matrix>,

   578             public random_access_iterator_base<packed_random_access_iterator_tag,

   579                                                const_iterator2, value_type> {

   580         public:

   581             typedef typename triangular_matrix::value_type value_type;

   582             typedef typename triangular_matrix::difference_type difference_type;

   583             typedef typename triangular_matrix::const_reference reference;

   584             typedef const typename triangular_matrix::pointer pointer;

   586             typedef const_iterator1 dual_iterator_type;

   587             typedef const_reverse_iterator1 dual_reverse_iterator_type;

   589             // Construction and destruction

   590             BOOST_UBLAS_INLINE

   591             const_iterator2 ():

   592                 container_const_reference<self_type> (), it1_ (), it2_ () {}

   593             BOOST_UBLAS_INLINE

   594             const_iterator2 (const self_type &m, size_type it1, size_type it2):

   595                 container_const_reference<self_type> (m), it1_ (it1), it2_ (it2) {}

   596             BOOST_UBLAS_INLINE

   597             const_iterator2 (const iterator2 &it):

   598                 container_const_reference<self_type> (it ()), it1_ (it.it1_), it2_ (it.it2_) {}

   600             // Arithmetic

   601             BOOST_UBLAS_INLINE

   602             const_iterator2 &operator ++ () {

   603                 ++ it2_;

   604                 return *this;

   605             }

   606             BOOST_UBLAS_INLINE

   607             const_iterator2 &operator -- () {

   608                 -- it2_;

   609                 return *this;

   610             }

   611             BOOST_UBLAS_INLINE

   612             const_iterator2 &operator += (difference_type n) {

   613                 it2_ += n;

   614                 return *this;

   615             }

   616             BOOST_UBLAS_INLINE

   617             const_iterator2 &operator -= (difference_type n) {

   618                 it2_ -= n;

   619                 return *this;

   620             }

   621             BOOST_UBLAS_INLINE

   622             difference_type operator - (const const_iterator2 &it) const {

   623                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

   624                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());

   625                 return it2_ - it.it2_;

   626             }

   628             // Dereference

   629             BOOST_UBLAS_INLINE

   630             const_reference operator * () const {

   631                 return (*this) () (it1_, it2_);

   632             }

   633             BOOST_UBLAS_INLINE

   634             const_reference operator [] (difference_type n) const {

   635                 return *(*this + n);

   636             }

   638 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION

   639             BOOST_UBLAS_INLINE

   640 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

   641             typename self_type::

   642 #endif

   643             const_iterator1 begin () const {

   644                 return (*this) ().find1 (1, 0, it2_);

   645             }

   646             BOOST_UBLAS_INLINE

   647 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

   648             typename self_type::

   649 #endif

   650             const_iterator1 end () const {

   651                 return (*this) ().find1 (1, (*this) ().size1 (), it2_);

   652             }

   653             BOOST_UBLAS_INLINE

   654 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

   655             typename self_type::

   656 #endif

   657             const_reverse_iterator1 rbegin () const {

   658                 return const_reverse_iterator1 (end ());

   659             }

   660             BOOST_UBLAS_INLINE

   661 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

   662             typename self_type::

   663 #endif

   664             const_reverse_iterator1 rend () const {

   665                 return const_reverse_iterator1 (begin ());

   666             }

   667 #endif

   669             // Indices

   670             BOOST_UBLAS_INLINE

   671             size_type index1 () const {

   672                 return it1_;

   673             }

   674             BOOST_UBLAS_INLINE

   675             size_type index2 () const {

   676                 return it2_;

   677             }

   679             // Assignment

   680             BOOST_UBLAS_INLINE

   681             const_iterator2 &operator = (const const_iterator2 &it) {

   682                 container_const_reference<self_type>::assign (&it ());

   683                 it1_ = it.it1_;

   684                 it2_ = it.it2_;

   685                 return *this;

   686             }

   688             // Comparison

   689             BOOST_UBLAS_INLINE

   690             bool operator == (const const_iterator2 &it) const {

   691                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

   692                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());

   693                 return it2_ == it.it2_;

   694             }

   695             BOOST_UBLAS_INLINE

   696             bool operator < (const const_iterator2 &it) const {

   697                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

   698                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());

   699                 return it2_ < it.it2_;

   700             }

   702         private:

   703             size_type it1_;

   704             size_type it2_;

   705         };

   706 #endif

   708         BOOST_UBLAS_INLINE

   709         const_iterator2 begin2 () const {

   710             return find2 (0, 0, 0);

   711         }

   712         BOOST_UBLAS_INLINE

   713         const_iterator2 end2 () const {

   714             return find2 (0, 0, size2_);

   715         }

   717 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR

   718         class iterator2:

   719             public container_reference<triangular_matrix>,

   720             public random_access_iterator_base<packed_random_access_iterator_tag,

   721                                                iterator2, value_type> {

   722         public:

   723             typedef typename triangular_matrix::value_type value_type;

   724             typedef typename triangular_matrix::difference_type difference_type;

   725             typedef typename triangular_matrix::reference reference;

   726             typedef typename triangular_matrix::pointer pointer;

   728             typedef iterator1 dual_iterator_type;

   729             typedef reverse_iterator1 dual_reverse_iterator_type;

   731             // Construction and destruction

   732             BOOST_UBLAS_INLINE

   733             iterator2 ():

   734                 container_reference<self_type> (), it1_ (), it2_ () {}

   735             BOOST_UBLAS_INLINE

   736             iterator2 (self_type &m, size_type it1, size_type it2):

   737                 container_reference<self_type> (m), it1_ (it1), it2_ (it2) {}

   739             // Arithmetic

   740             BOOST_UBLAS_INLINE

   741             iterator2 &operator ++ () {

   742                 ++ it2_;

   743                 return *this;

   744             }

   745             BOOST_UBLAS_INLINE

   746             iterator2 &operator -- () {

   747                 -- it2_;

   748                 return *this;

   749             }

   750             BOOST_UBLAS_INLINE

   751             iterator2 &operator += (difference_type n) {

   752                 it2_ += n;

   753                 return *this;

   754             }

   755             BOOST_UBLAS_INLINE

   756             iterator2 &operator -= (difference_type n) {

   757                 it2_ -= n;

   758                 return *this;

   759             }

   760             BOOST_UBLAS_INLINE

   761             difference_type operator - (const iterator2 &it) const {

   762                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

   763                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());

   764                 return it2_ - it.it2_;

   765             }

   767             // Dereference

   768             BOOST_UBLAS_INLINE

   769             reference operator * () const {

   770                 return (*this) () (it1_, it2_);

   771             }

   772             BOOST_UBLAS_INLINE

   773             reference operator [] (difference_type n) const {

   774                 return *(*this + n);

   775             }

   777 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION

   778             BOOST_UBLAS_INLINE

   779 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

   780             typename self_type::

   781 #endif

   782             iterator1 begin () const {

   783                 return (*this) ().find1 (1, 0, it2_);

   784             }

   785             BOOST_UBLAS_INLINE

   786 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

   787             typename self_type::

   788 #endif

   789             iterator1 end () const {

   790                 return (*this) ().find1 (1, (*this) ().size1 (), it2_);

   791             }

   792             BOOST_UBLAS_INLINE

   793 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

   794             typename self_type::

   795 #endif

   796             reverse_iterator1 rbegin () const {

   797                 return reverse_iterator1 (end ());

   798             }

   799             BOOST_UBLAS_INLINE

   800 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

   801             typename self_type::

   802 #endif

   803             reverse_iterator1 rend () const {

   804                 return reverse_iterator1 (begin ());

   805             }

   806 #endif

   808             // Indices

   809             BOOST_UBLAS_INLINE

   810             size_type index1 () const {

   811                 return it1_;

   812             }

   813             BOOST_UBLAS_INLINE

   814             size_type index2 () const {

   815                 return it2_;

   816             }

   818             // Assignment

   819             BOOST_UBLAS_INLINE

   820             iterator2 &operator = (const iterator2 &it) {

   821                 container_reference<self_type>::assign (&it ());

   822                 it1_ = it.it1_;

   823                 it2_ = it.it2_;

   824                 return *this;

   825             }

   827             // Comparison

   828             BOOST_UBLAS_INLINE

   829             bool operator == (const iterator2 &it) const {

   830                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

   831                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());

   832                 return it2_ == it.it2_;

   833             }

   834             BOOST_UBLAS_INLINE

   835             bool operator < (const iterator2 &it) const {

   836                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

   837                 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());

   838                 return it2_ < it.it2_;

   839             }

   841         private:

   842             size_type it1_;

   843             size_type it2_;

   845             friend class const_iterator2;

   846         };

   847 #endif

   849         BOOST_UBLAS_INLINE

   850         iterator2 begin2 () {

   851             return find2 (0, 0, 0);

   852         }

   853         BOOST_UBLAS_INLINE

   854         iterator2 end2 () {

   855             return find2 (0, 0, size2_);

   856         }

   858         // Reverse iterators

   860         BOOST_UBLAS_INLINE

   861         const_reverse_iterator1 rbegin1 () const {

   862             return const_reverse_iterator1 (end1 ());

   863         }

   864         BOOST_UBLAS_INLINE

   865         const_reverse_iterator1 rend1 () const {

   866             return const_reverse_iterator1 (begin1 ());

   867         }

   869         BOOST_UBLAS_INLINE

   870         reverse_iterator1 rbegin1 () {

   871             return reverse_iterator1 (end1 ());

   872         }

   873         BOOST_UBLAS_INLINE

   874         reverse_iterator1 rend1 () {

   875             return reverse_iterator1 (begin1 ());

   876         }

   878         BOOST_UBLAS_INLINE

   879         const_reverse_iterator2 rbegin2 () const {

   880             return const_reverse_iterator2 (end2 ());

   881         }

   882         BOOST_UBLAS_INLINE

   883         const_reverse_iterator2 rend2 () const {

   884             return const_reverse_iterator2 (begin2 ());

   885         }

   887         BOOST_UBLAS_INLINE

   888         reverse_iterator2 rbegin2 () {

   889             return reverse_iterator2 (end2 ());

   890         }

   891         BOOST_UBLAS_INLINE

   892         reverse_iterator2 rend2 () {

   893             return reverse_iterator2 (begin2 ());

   894         }

   896     private:

   897         size_type size1_;

   898         size_type size2_;

   899         array_type data_;

   900         static const value_type zero_;

   901         static const value_type one_;

   902     };

   904     template<class T, class TRI, class L, class A>

   905     const typename triangular_matrix<T, TRI, L, A>::value_type triangular_matrix<T, TRI, L, A>::zero_ = value_type/*zero*/();

   906     template<class T, class TRI, class L, class A>

   907     const typename triangular_matrix<T, TRI, L, A>::value_type triangular_matrix<T, TRI, L, A>::one_ (1);

   910     // Triangular matrix adaptor class

   911     template<class M, class TRI>

   912     class triangular_adaptor:

   913         public matrix_expression<triangular_adaptor<M, TRI> > {

   915         typedef triangular_adaptor<M, TRI> self_type;

   917     public:

   918 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS

   919         using matrix_expression<self_type>::operator ();

   920 #endif

   921         typedef const M const_matrix_type;

   922         typedef M matrix_type;

   923         typedef TRI triangular_type;

   924         typedef typename M::size_type size_type;

   925         typedef typename M::difference_type difference_type;

   926         typedef typename M::value_type value_type;

   927         typedef typename M::const_reference const_reference;

   928         typedef typename boost::mpl::if_<boost::is_const<M>,

   929                                           typename M::const_reference,

   930                                           typename M::reference>::type reference;

   931         typedef typename boost::mpl::if_<boost::is_const<M>,

   932                                           typename M::const_closure_type,

   933                                           typename M::closure_type>::type matrix_closure_type;

   934         typedef const self_type const_closure_type;

   935         typedef self_type closure_type;

   936         // Replaced by _temporary_traits to avoid type requirements on M

   937         //typedef typename M::vector_temporary_type vector_temporary_type;

   938         //typedef typename M::matrix_temporary_type matrix_temporary_type;

   939         typedef typename storage_restrict_traits<typename M::storage_category,

   940                                                  packed_proxy_tag>::storage_category storage_category;

   941         typedef typename M::orientation_category orientation_category;

   943         // Construction and destruction

   944         BOOST_UBLAS_INLINE

   945         triangular_adaptor (matrix_type &data):

   946             matrix_expression<self_type> (),

   947             data_ (data) {}

   948         BOOST_UBLAS_INLINE

   949         triangular_adaptor (const triangular_adaptor &m):

   950             matrix_expression<self_type> (),

   951             data_ (m.data_) {}

   953         // Accessors

   954         BOOST_UBLAS_INLINE

   955         size_type size1 () const {

   956             return data_.size1 ();

   957         }

   958         BOOST_UBLAS_INLINE

   959         size_type size2 () const {

   960             return data_.size2 ();

   961         }

   963         // Storage accessors

   964         BOOST_UBLAS_INLINE

   965         const matrix_closure_type &data () const {

   966             return data_;

   967         }

   968         BOOST_UBLAS_INLINE

   969         matrix_closure_type &data () {

   970             return data_;

   971         }

   973         // Element access

   974 #ifndef BOOST_UBLAS_PROXY_CONST_MEMBER

   975         BOOST_UBLAS_INLINE

   976         const_reference operator () (size_type i, size_type j) const {

   977             BOOST_UBLAS_CHECK (i < size1 (), bad_index ());

   978             BOOST_UBLAS_CHECK (j < size2 (), bad_index ());

   979             if (triangular_type::other (i, j))

   980                 return data () (i, j);

   981             else if (triangular_type::one (i, j))

   982                 return one_;

   983             else

   984                 return zero_;

   985         }

   986         BOOST_UBLAS_INLINE

   987         reference operator () (size_type i, size_type j) {

   988             BOOST_UBLAS_CHECK (i < size1 (), bad_index ());

   989             BOOST_UBLAS_CHECK (j < size2 (), bad_index ());

   990             if (triangular_type::other (i, j))

   991                 return data () (i, j);

   992             else if (triangular_type::one (i, j)) {

   993                 bad_index ().raise ();

   994                 // arbitary return value

   995                 return const_cast<reference>(one_);

   996             } else {

   997                 bad_index ().raise ();

   998                 // arbitary return value

   999                 return const_cast<reference>(zero_);

  1000             }

  1001         }

  1002 #else

  1003         BOOST_UBLAS_INLINE

  1004         reference operator () (size_type i, size_type j) const {

  1005             BOOST_UBLAS_CHECK (i < size1 (), bad_index ());

  1006             BOOST_UBLAS_CHECK (j < size2 (), bad_index ());

  1007             if (triangular_type::other (i, j))

  1008                 return data () (i, j);

  1009             else if (triangular_type::one (i, j)) {

  1010                 bad_index ().raise ();

  1011                 // arbitary return value

  1012                 return const_cast<reference>(one_);

  1013             } else {

  1014                 bad_index ().raise ();

  1015                 // arbitary return value

  1016                 return const_cast<reference>(zero_);

  1017             }

  1018         }

  1019 #endif

  1021         // Assignment

  1022         BOOST_UBLAS_INLINE

  1023         triangular_adaptor &operator = (const triangular_adaptor &m) {

  1024             matrix_assign<scalar_assign> (*this, m);

  1025             return *this;

  1026         }

  1027         BOOST_UBLAS_INLINE

  1028         triangular_adaptor &assign_temporary (triangular_adaptor &m) {

  1029             *this = m;

  1030             return *this;

  1031         }

  1032         template<class AE>

  1033         BOOST_UBLAS_INLINE

  1034         triangular_adaptor &operator = (const matrix_expression<AE> &ae) {

  1035             matrix_assign<scalar_assign> (*this, matrix<value_type> (ae));

  1036             return *this;

  1037         }

  1038         template<class AE>

  1039         BOOST_UBLAS_INLINE

  1040         triangular_adaptor &assign (const matrix_expression<AE> &ae) {

  1041             matrix_assign<scalar_assign> (*this, ae);

  1042             return *this;

  1043         }

  1044         template<class AE>

  1045         BOOST_UBLAS_INLINE

  1046         triangular_adaptor& operator += (const matrix_expression<AE> &ae) {

  1047             matrix_assign<scalar_assign> (*this, matrix<value_type> (*this + ae));

  1048             return *this;

  1049         }

  1050         template<class AE>

  1051         BOOST_UBLAS_INLINE

  1052         triangular_adaptor &plus_assign (const matrix_expression<AE> &ae) {

  1053             matrix_assign<scalar_plus_assign> (*this, ae);

  1054             return *this;

  1055         }

  1056         template<class AE>

  1057         BOOST_UBLAS_INLINE

  1058         triangular_adaptor& operator -= (const matrix_expression<AE> &ae) {

  1059             matrix_assign<scalar_assign> (*this, matrix<value_type> (*this - ae));

  1060             return *this;

  1061         }

  1062         template<class AE>

  1063         BOOST_UBLAS_INLINE

  1064         triangular_adaptor &minus_assign (const matrix_expression<AE> &ae) {

  1065             matrix_assign<scalar_minus_assign> (*this, ae);

  1066             return *this;

  1067         }

  1068         template<class AT>

  1069         BOOST_UBLAS_INLINE

  1070         triangular_adaptor& operator *= (const AT &at) {

  1071             matrix_assign_scalar<scalar_multiplies_assign> (*this, at);

  1072             return *this;

  1073         }

  1074         template<class AT>

  1075         BOOST_UBLAS_INLINE

  1076         triangular_adaptor& operator /= (const AT &at) {

  1077             matrix_assign_scalar<scalar_divides_assign> (*this, at);

  1078             return *this;

  1079         }

  1081         // Closure comparison

  1082         BOOST_UBLAS_INLINE

  1083         bool same_closure (const triangular_adaptor &ta) const {

  1084             return (*this).data ().same_closure (ta.data ());

  1085         }

  1087         // Swapping

  1088         BOOST_UBLAS_INLINE

  1089         void swap (triangular_adaptor &m) {

  1090             if (this != &m)

  1091                 matrix_swap<scalar_swap> (*this, m);

  1092         }

  1093         BOOST_UBLAS_INLINE

  1094         friend void swap (triangular_adaptor &m1, triangular_adaptor &m2) {

  1095             m1.swap (m2);

  1096         }

  1098         // Iterator types

  1099    private:

  1100         typedef typename M::const_iterator1 const_subiterator1_type;

  1101         typedef typename boost::mpl::if_<boost::is_const<M>,

  1102                                           typename M::const_iterator1,

  1103                                           typename M::iterator1>::type subiterator1_type;

  1104         typedef typename M::const_iterator2 const_subiterator2_type;

  1105         typedef typename boost::mpl::if_<boost::is_const<M>,

  1106                                           typename M::const_iterator2,

  1107                                           typename M::iterator2>::type subiterator2_type;

  1109     public:

  1110 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR

  1111         typedef indexed_iterator1<self_type, packed_random_access_iterator_tag> iterator1;

  1112         typedef indexed_iterator2<self_type, packed_random_access_iterator_tag> iterator2;

  1113         typedef indexed_const_iterator1<self_type, packed_random_access_iterator_tag> const_iterator1;

  1114         typedef indexed_const_iterator2<self_type, packed_random_access_iterator_tag> const_iterator2;

  1115 #else

  1116         class const_iterator1;

  1117         class iterator1;

  1118         class const_iterator2;

  1119         class iterator2;

  1120 #endif

  1121         typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;

  1122         typedef reverse_iterator_base1<iterator1> reverse_iterator1;

  1123         typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;

  1124         typedef reverse_iterator_base2<iterator2> reverse_iterator2;

  1126         // Element lookup

  1127         BOOST_UBLAS_INLINE

  1128         const_iterator1 find1 (int rank, size_type i, size_type j) const {

  1129             if (rank == 1)

  1130                 i = triangular_type::restrict1 (i, j);

  1131             return const_iterator1 (*this, data ().find1 (rank, i, j));

  1132         }

  1133         BOOST_UBLAS_INLINE

  1134         iterator1 find1 (int rank, size_type i, size_type j) {

  1135             if (rank == 1)

  1136                 i = triangular_type::mutable_restrict1 (i, j);

  1137             return iterator1 (*this, data ().find1 (rank, i, j));

  1138         }

  1139         BOOST_UBLAS_INLINE

  1140         const_iterator2 find2 (int rank, size_type i, size_type j) const {

  1141             if (rank == 1)

  1142                 j = triangular_type::restrict2 (i, j);

  1143             return const_iterator2 (*this, data ().find2 (rank, i, j));

  1144         }

  1145         BOOST_UBLAS_INLINE

  1146         iterator2 find2 (int rank, size_type i, size_type j) {

  1147             if (rank == 1)

  1148                 j = triangular_type::mutable_restrict2 (i, j);

  1149             return iterator2 (*this, data ().find2 (rank, i, j));

  1150         }

  1152         // Iterators simply are indices.

  1154 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR

  1155         class const_iterator1:

  1156             public container_const_reference<triangular_adaptor>,

  1157             public random_access_iterator_base<typename iterator_restrict_traits<

  1158                                                    typename const_subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,

  1159                                                const_iterator1, value_type> {

  1160         public:

  1161             typedef typename const_subiterator1_type::value_type value_type;

  1162             typedef typename const_subiterator1_type::difference_type difference_type;

  1163             typedef typename const_subiterator1_type::reference reference;

  1164             typedef typename const_subiterator1_type::pointer pointer;

  1166             typedef const_iterator2 dual_iterator_type;

  1167             typedef const_reverse_iterator2 dual_reverse_iterator_type;

  1169             // Construction and destruction

  1170             BOOST_UBLAS_INLINE

  1171             const_iterator1 ():

  1172                 container_const_reference<self_type> (), it1_ () {}

  1173             BOOST_UBLAS_INLINE

  1174             const_iterator1 (const self_type &m, const const_subiterator1_type &it1):

  1175                 container_const_reference<self_type> (m), it1_ (it1) {}

  1176             BOOST_UBLAS_INLINE

  1177             const_iterator1 (const iterator1 &it):

  1178                 container_const_reference<self_type> (it ()), it1_ (it.it1_) {}

  1180             // Arithmetic

  1181             BOOST_UBLAS_INLINE

  1182             const_iterator1 &operator ++ () {

  1183                 ++ it1_;

  1184                 return *this;

  1185             }

  1186             BOOST_UBLAS_INLINE

  1187             const_iterator1 &operator -- () {

  1188                 -- it1_;

  1189                 return *this;

  1190             }

  1191             BOOST_UBLAS_INLINE

  1192             const_iterator1 &operator += (difference_type n) {

  1193                 it1_ += n;

  1194                 return *this;

  1195             }

  1196             BOOST_UBLAS_INLINE

  1197             const_iterator1 &operator -= (difference_type n) {

  1198                 it1_ -= n;

  1199                 return *this;

  1200             }

  1201             BOOST_UBLAS_INLINE

  1202             difference_type operator - (const const_iterator1 &it) const {

  1203                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

  1204                 return it1_ - it.it1_;

  1205             }

  1207             // Dereference

  1208             BOOST_UBLAS_INLINE

  1209             const_reference operator * () const {

  1210                 size_type i = index1 ();

  1211                 size_type j = index2 ();

  1212                 BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());

  1213                 BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());

  1214                 if (triangular_type::other (i, j))

  1215                     return *it1_;

  1216                 else

  1217                     return (*this) () (i, j);

  1218             }

  1219             BOOST_UBLAS_INLINE

  1220             const_reference operator [] (difference_type n) const {

  1221                 return *(*this + n);

  1222             }

  1224 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION

  1225             BOOST_UBLAS_INLINE

  1226 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

  1227             typename self_type::

  1228 #endif

  1229             const_iterator2 begin () const {

  1230                 return (*this) ().find2 (1, index1 (), 0);

  1231             }

  1232             BOOST_UBLAS_INLINE

  1233 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

  1234             typename self_type::

  1235 #endif

  1236             const_iterator2 end () const {

  1237                 return (*this) ().find2 (1, index1 (), (*this) ().size2 ());

  1238             }

  1239             BOOST_UBLAS_INLINE

  1240 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

  1241             typename self_type::

  1242 #endif

  1243             const_reverse_iterator2 rbegin () const {

  1244                 return const_reverse_iterator2 (end ());

  1245             }

  1246             BOOST_UBLAS_INLINE

  1247 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

  1248             typename self_type::

  1249 #endif

  1250             const_reverse_iterator2 rend () const {

  1251                 return const_reverse_iterator2 (begin ());

  1252             }

  1253 #endif

  1255             // Indices

  1256             BOOST_UBLAS_INLINE

  1257             size_type index1 () const {

  1258                 return it1_.index1 ();

  1259             }

  1260             BOOST_UBLAS_INLINE

  1261             size_type index2 () const {

  1262                 return it1_.index2 ();

  1263             }

  1265             // Assignment

  1266             BOOST_UBLAS_INLINE

  1267             const_iterator1 &operator = (const const_iterator1 &it) {

  1268                 container_const_reference<self_type>::assign (&it ());

  1269                 it1_ = it.it1_;

  1270                 return *this;

  1271             }

  1273             // Comparison

  1274             BOOST_UBLAS_INLINE

  1275             bool operator == (const const_iterator1 &it) const {

  1276                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

  1277                 return it1_ == it.it1_;

  1278             }

  1279             BOOST_UBLAS_INLINE

  1280             bool operator < (const const_iterator1 &it) const {

  1281                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

  1282                 return it1_ < it.it1_;

  1283             }

  1285         private:

  1286             const_subiterator1_type it1_;

  1287         };

  1288 #endif

  1290         BOOST_UBLAS_INLINE

  1291         const_iterator1 begin1 () const {

  1292             return find1 (0, 0, 0);

  1293         }

  1294         BOOST_UBLAS_INLINE

  1295         const_iterator1 end1 () const {

  1296             return find1 (0, size1 (), 0);

  1297         }

  1299 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR

  1300         class iterator1:

  1301             public container_reference<triangular_adaptor>,

  1302             public random_access_iterator_base<typename iterator_restrict_traits<

  1303                                                    typename subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,

  1304                                                iterator1, value_type> {

  1305         public:

  1306             typedef typename subiterator1_type::value_type value_type;

  1307             typedef typename subiterator1_type::difference_type difference_type;

  1308             typedef typename subiterator1_type::reference reference;

  1309             typedef typename subiterator1_type::pointer pointer;

  1311             typedef iterator2 dual_iterator_type;

  1312             typedef reverse_iterator2 dual_reverse_iterator_type;

  1314             // Construction and destruction

  1315             BOOST_UBLAS_INLINE

  1316             iterator1 ():

  1317                 container_reference<self_type> (), it1_ () {}

  1318             BOOST_UBLAS_INLINE

  1319             iterator1 (self_type &m, const subiterator1_type &it1):

  1320                 container_reference<self_type> (m), it1_ (it1) {}

  1322             // Arithmetic

  1323             BOOST_UBLAS_INLINE

  1324             iterator1 &operator ++ () {

  1325                 ++ it1_;

  1326                 return *this;

  1327             }

  1328             BOOST_UBLAS_INLINE

  1329             iterator1 &operator -- () {

  1330                 -- it1_;

  1331                 return *this;

  1332             }

  1333             BOOST_UBLAS_INLINE

  1334             iterator1 &operator += (difference_type n) {

  1335                 it1_ += n;

  1336                 return *this;

  1337             }

  1338             BOOST_UBLAS_INLINE

  1339             iterator1 &operator -= (difference_type n) {

  1340                 it1_ -= n;

  1341                 return *this;

  1342             }

  1343             BOOST_UBLAS_INLINE

  1344             difference_type operator - (const iterator1 &it) const {

  1345                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

  1346                 return it1_ - it.it1_;

  1347             }

  1349             // Dereference

  1350             BOOST_UBLAS_INLINE

  1351             reference operator * () const {

  1352                 size_type i = index1 ();

  1353                 size_type j = index2 ();

  1354                 BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());

  1355                 BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());

  1356                 if (triangular_type::other (i, j))

  1357                     return *it1_;

  1358                 else

  1359                     return (*this) () (i, j);

  1360             }

  1361             BOOST_UBLAS_INLINE

  1362             reference operator [] (difference_type n) const {

  1363                 return *(*this + n);

  1364             }

  1366 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION

  1367             BOOST_UBLAS_INLINE

  1368 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

  1369             typename self_type::

  1370 #endif

  1371             iterator2 begin () const {

  1372                 return (*this) ().find2 (1, index1 (), 0);

  1373             }

  1374             BOOST_UBLAS_INLINE

  1375 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

  1376             typename self_type::

  1377 #endif

  1378             iterator2 end () const {

  1379                 return (*this) ().find2 (1, index1 (), (*this) ().size2 ());

  1380             }

  1381             BOOST_UBLAS_INLINE

  1382 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

  1383             typename self_type::

  1384 #endif

  1385             reverse_iterator2 rbegin () const {

  1386                 return reverse_iterator2 (end ());

  1387             }

  1388             BOOST_UBLAS_INLINE

  1389 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

  1390             typename self_type::

  1391 #endif

  1392             reverse_iterator2 rend () const {

  1393                 return reverse_iterator2 (begin ());

  1394             }

  1395 #endif

  1397             // Indices

  1398             BOOST_UBLAS_INLINE

  1399             size_type index1 () const {

  1400                 return it1_.index1 ();

  1401             }

  1402             BOOST_UBLAS_INLINE

  1403             size_type index2 () const {

  1404                 return it1_.index2 ();

  1405             }

  1407             // Assignment

  1408             BOOST_UBLAS_INLINE

  1409             iterator1 &operator = (const iterator1 &it) {

  1410                 container_reference<self_type>::assign (&it ());

  1411                 it1_ = it.it1_;

  1412                 return *this;

  1413             }

  1415             // Comparison

  1416             BOOST_UBLAS_INLINE

  1417             bool operator == (const iterator1 &it) const {

  1418                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

  1419                 return it1_ == it.it1_;

  1420             }

  1421             BOOST_UBLAS_INLINE

  1422             bool operator < (const iterator1 &it) const {

  1423                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

  1424                 return it1_ < it.it1_;

  1425             }

  1427         private:

  1428             subiterator1_type it1_;

  1430             friend class const_iterator1;

  1431         };

  1432 #endif

  1434         BOOST_UBLAS_INLINE

  1435         iterator1 begin1 () {

  1436             return find1 (0, 0, 0);

  1437         }

  1438         BOOST_UBLAS_INLINE

  1439         iterator1 end1 () {

  1440             return find1 (0, size1 (), 0);

  1441         }

  1443 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR

  1444         class const_iterator2:

  1445             public container_const_reference<triangular_adaptor>,

  1446             public random_access_iterator_base<typename iterator_restrict_traits<

  1447                                                    typename const_subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,

  1448                                                const_iterator2, value_type> {

  1449         public:

  1450             typedef typename const_subiterator2_type::value_type value_type;

  1451             typedef typename const_subiterator2_type::difference_type difference_type;

  1452             typedef typename const_subiterator2_type::reference reference;

  1453             typedef typename const_subiterator2_type::pointer pointer;

  1455             typedef const_iterator1 dual_iterator_type;

  1456             typedef const_reverse_iterator1 dual_reverse_iterator_type;

  1458             // Construction and destruction

  1459             BOOST_UBLAS_INLINE

  1460             const_iterator2 ():

  1461                 container_const_reference<self_type> (), it2_ () {}

  1462             BOOST_UBLAS_INLINE

  1463             const_iterator2 (const self_type &m, const const_subiterator2_type &it2):

  1464                 container_const_reference<self_type> (m), it2_ (it2) {}

  1465             BOOST_UBLAS_INLINE

  1466             const_iterator2 (const iterator2 &it):

  1467                 container_const_reference<self_type> (it ()), it2_ (it.it2_) {}

  1469             // Arithmetic

  1470             BOOST_UBLAS_INLINE

  1471             const_iterator2 &operator ++ () {

  1472                 ++ it2_;

  1473                 return *this;

  1474             }

  1475             BOOST_UBLAS_INLINE

  1476             const_iterator2 &operator -- () {

  1477                 -- it2_;

  1478                 return *this;

  1479             }

  1480             BOOST_UBLAS_INLINE

  1481             const_iterator2 &operator += (difference_type n) {

  1482                 it2_ += n;

  1483                 return *this;

  1484             }

  1485             BOOST_UBLAS_INLINE

  1486             const_iterator2 &operator -= (difference_type n) {

  1487                 it2_ -= n;

  1488                 return *this;

  1489             }

  1490             BOOST_UBLAS_INLINE

  1491             difference_type operator - (const const_iterator2 &it) const {

  1492                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

  1493                 return it2_ - it.it2_;

  1494             }

  1496             // Dereference

  1497             BOOST_UBLAS_INLINE

  1498             const_reference operator * () const {

  1499                 size_type i = index1 ();

  1500                 size_type j = index2 ();

  1501                 BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());

  1502                 BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());

  1503                 if (triangular_type::other (i, j))

  1504                     return *it2_;

  1505                 else

  1506                     return (*this) () (i, j);

  1507             }

  1508             BOOST_UBLAS_INLINE

  1509             const_reference operator [] (difference_type n) const {

  1510                 return *(*this + n);

  1511             }

  1513 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION

  1514             BOOST_UBLAS_INLINE

  1515 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

  1516             typename self_type::

  1517 #endif

  1518             const_iterator1 begin () const {

  1519                 return (*this) ().find1 (1, 0, index2 ());

  1520             }

  1521             BOOST_UBLAS_INLINE

  1522 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

  1523             typename self_type::

  1524 #endif

  1525             const_iterator1 end () const {

  1526                 return (*this) ().find1 (1, (*this) ().size1 (), index2 ());

  1527             }

  1528             BOOST_UBLAS_INLINE

  1529 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

  1530             typename self_type::

  1531 #endif

  1532             const_reverse_iterator1 rbegin () const {

  1533                 return const_reverse_iterator1 (end ());

  1534             }

  1535             BOOST_UBLAS_INLINE

  1536 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

  1537             typename self_type::

  1538 #endif

  1539             const_reverse_iterator1 rend () const {

  1540                 return const_reverse_iterator1 (begin ());

  1541             }

  1542 #endif

  1544             // Indices

  1545             BOOST_UBLAS_INLINE

  1546             size_type index1 () const {

  1547                 return it2_.index1 ();

  1548             }

  1549             BOOST_UBLAS_INLINE

  1550             size_type index2 () const {

  1551                 return it2_.index2 ();

  1552             }

  1554             // Assignment

  1555             BOOST_UBLAS_INLINE

  1556             const_iterator2 &operator = (const const_iterator2 &it) {

  1557                 container_const_reference<self_type>::assign (&it ());

  1558                 it2_ = it.it2_;

  1559                 return *this;

  1560             }

  1562             // Comparison

  1563             BOOST_UBLAS_INLINE

  1564             bool operator == (const const_iterator2 &it) const {

  1565                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

  1566                 return it2_ == it.it2_;

  1567             }

  1568             BOOST_UBLAS_INLINE

  1569             bool operator < (const const_iterator2 &it) const {

  1570                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

  1571                 return it2_ < it.it2_;

  1572             }

  1574         private:

  1575             const_subiterator2_type it2_;

  1576         };

  1577 #endif

  1579         BOOST_UBLAS_INLINE

  1580         const_iterator2 begin2 () const {

  1581             return find2 (0, 0, 0);

  1582         }

  1583         BOOST_UBLAS_INLINE

  1584         const_iterator2 end2 () const {

  1585             return find2 (0, 0, size2 ());

  1586         }

  1588 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR

  1589         class iterator2:

  1590             public container_reference<triangular_adaptor>,

  1591             public random_access_iterator_base<typename iterator_restrict_traits<

  1592                                                    typename subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,

  1593                                                iterator2, value_type> {

  1594         public:

  1595             typedef typename subiterator2_type::value_type value_type;

  1596             typedef typename subiterator2_type::difference_type difference_type;

  1597             typedef typename subiterator2_type::reference reference;

  1598             typedef typename subiterator2_type::pointer pointer;

  1600             typedef iterator1 dual_iterator_type;

  1601             typedef reverse_iterator1 dual_reverse_iterator_type;

  1603             // Construction and destruction

  1604             BOOST_UBLAS_INLINE

  1605             iterator2 ():

  1606                 container_reference<self_type> (), it2_ () {}

  1607             BOOST_UBLAS_INLINE

  1608             iterator2 (self_type &m, const subiterator2_type &it2):

  1609                 container_reference<self_type> (m), it2_ (it2) {}

  1611             // Arithmetic

  1612             BOOST_UBLAS_INLINE

  1613             iterator2 &operator ++ () {

  1614                 ++ it2_;

  1615                 return *this;

  1616             }

  1617             BOOST_UBLAS_INLINE

  1618             iterator2 &operator -- () {

  1619                 -- it2_;

  1620                 return *this;

  1621             }

  1622             BOOST_UBLAS_INLINE

  1623             iterator2 &operator += (difference_type n) {

  1624                 it2_ += n;

  1625                 return *this;

  1626             }

  1627             BOOST_UBLAS_INLINE

  1628             iterator2 &operator -= (difference_type n) {

  1629                 it2_ -= n;

  1630                 return *this;

  1631             }

  1632             BOOST_UBLAS_INLINE

  1633             difference_type operator - (const iterator2 &it) const {

  1634                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

  1635                 return it2_ - it.it2_;

  1636             }

  1638             // Dereference

  1639             BOOST_UBLAS_INLINE

  1640             reference operator * () const {

  1641                 size_type i = index1 ();

  1642                 size_type j = index2 ();

  1643                 BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());

  1644                 BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());

  1645                 if (triangular_type::other (i, j))

  1646                     return *it2_;

  1647                 else

  1648                     return (*this) () (i, j);

  1649             }

  1650             BOOST_UBLAS_INLINE

  1651             reference operator [] (difference_type n) const {

  1652                 return *(*this + n);

  1653             }

  1655 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION

  1656             BOOST_UBLAS_INLINE

  1657 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

  1658             typename self_type::

  1659 #endif

  1660             iterator1 begin () const {

  1661                 return (*this) ().find1 (1, 0, index2 ());

  1662             }

  1663             BOOST_UBLAS_INLINE

  1664 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

  1665             typename self_type::

  1666 #endif

  1667             iterator1 end () const {

  1668                 return (*this) ().find1 (1, (*this) ().size1 (), index2 ());

  1669             }

  1670             BOOST_UBLAS_INLINE

  1671 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

  1672             typename self_type::

  1673 #endif

  1674             reverse_iterator1 rbegin () const {

  1675                 return reverse_iterator1 (end ());

  1676             }

  1677             BOOST_UBLAS_INLINE

  1678 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION

  1679             typename self_type::

  1680 #endif

  1681             reverse_iterator1 rend () const {

  1682                 return reverse_iterator1 (begin ());

  1683             }

  1684 #endif

  1686             // Indices

  1687             BOOST_UBLAS_INLINE

  1688             size_type index1 () const {

  1689                 return it2_.index1 ();

  1690             }

  1691             BOOST_UBLAS_INLINE

  1692             size_type index2 () const {

  1693                 return it2_.index2 ();

  1694             }

  1696             // Assignment

  1697             BOOST_UBLAS_INLINE

  1698             iterator2 &operator = (const iterator2 &it) {

  1699                 container_reference<self_type>::assign (&it ());

  1700                 it2_ = it.it2_;

  1701                 return *this;

  1702             }

  1704             // Comparison

  1705             BOOST_UBLAS_INLINE

  1706             bool operator == (const iterator2 &it) const {

  1707                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

  1708                 return it2_ == it.it2_;

  1709             }

  1710             BOOST_UBLAS_INLINE

  1711             bool operator < (const iterator2 &it) const {

  1712                 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());

  1713                 return it2_ < it.it2_;

  1714             }

  1716         private:

  1717             subiterator2_type it2_;

  1719             friend class const_iterator2;

  1720         };

  1721 #endif

  1723         BOOST_UBLAS_INLINE

  1724         iterator2 begin2 () {

  1725             return find2 (0, 0, 0);

  1726         }

  1727         BOOST_UBLAS_INLINE

  1728         iterator2 end2 () {

  1729             return find2 (0, 0, size2 ());

  1730         }

  1732         // Reverse iterators

  1734         BOOST_UBLAS_INLINE

  1735         const_reverse_iterator1 rbegin1 () const {

  1736             return const_reverse_iterator1 (end1 ());

  1737         }

  1738         BOOST_UBLAS_INLINE

  1739         const_reverse_iterator1 rend1 () const {

  1740             return const_reverse_iterator1 (begin1 ());

  1741         }

  1743         BOOST_UBLAS_INLINE

  1744         reverse_iterator1 rbegin1 () {

  1745             return reverse_iterator1 (end1 ());

  1746         }

  1747         BOOST_UBLAS_INLINE

  1748         reverse_iterator1 rend1 () {

  1749             return reverse_iterator1 (begin1 ());

  1750         }

  1752         BOOST_UBLAS_INLINE

  1753         const_reverse_iterator2 rbegin2 () const {

  1754             return const_reverse_iterator2 (end2 ());

  1755         }

  1756         BOOST_UBLAS_INLINE

  1757         const_reverse_iterator2 rend2 () const {

  1758             return const_reverse_iterator2 (begin2 ());

  1759         }

  1761         BOOST_UBLAS_INLINE

  1762         reverse_iterator2 rbegin2 () {

  1763             return reverse_iterator2 (end2 ());

  1764         }

  1765         BOOST_UBLAS_INLINE

  1766         reverse_iterator2 rend2 () {

  1767             return reverse_iterator2 (begin2 ());

  1768         }

  1770     private:

  1771         matrix_closure_type data_;

  1772         static const value_type zero_;

  1773         static const value_type one_;

  1774     };

  1776     template<class M, class TRI>

  1777     const typename triangular_adaptor<M, TRI>::value_type triangular_adaptor<M, TRI>::zero_ = value_type/*zero*/();

  1778     template<class M, class TRI>

  1779     const typename triangular_adaptor<M, TRI>::value_type triangular_adaptor<M, TRI>::one_ (1);

  1781     template <class M, class TRI>

  1782     struct vector_temporary_traits< triangular_adaptor<M, TRI> >

  1783     : vector_temporary_traits< typename boost::remove_const<M>::type > {} ;

  1784     template <class M, class TRI>

  1785     struct vector_temporary_traits< const triangular_adaptor<M, TRI> >

  1786     : vector_temporary_traits< typename boost::remove_const<M>::type > {} ;

  1788     template <class M, class TRI>

  1789     struct matrix_temporary_traits< triangular_adaptor<M, TRI> >

  1790     : matrix_temporary_traits< typename boost::remove_const<M>::type > {};

  1791     template <class M, class TRI>

  1792     struct matrix_temporary_traits< const triangular_adaptor<M, TRI> >

  1793     : matrix_temporary_traits< typename boost::remove_const<M>::type > {};

  1796     template<class E1, class E2>

  1797     struct matrix_vector_solve_traits {

  1798         typedef typename promote_traits<typename E1::value_type, typename E2::value_type>::promote_type promote_type;

  1799         typedef vector<promote_type> result_type;

  1800     };

  1802     // Operations:

  1803     //  n * (n - 1) / 2 + n = n * (n + 1) / 2 multiplications,

  1804     //  n * (n - 1) / 2 additions

  1806     // Dense (proxy) case

  1807     template<class E1, class E2>

  1808     BOOST_UBLAS_INLINE

  1809     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,

  1810                         lower_tag, column_major_tag, dense_proxy_tag) {

  1811         typedef typename E2::size_type size_type;

  1812         typedef typename E2::difference_type difference_type;

  1813         typedef typename E2::value_type value_type;

  1815         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());

  1816         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());

  1817         size_type size = e2 ().size ();

  1818         for (size_type n = 0; n < size; ++ n) {

  1819 #ifndef BOOST_UBLAS_SINGULAR_CHECK

  1820             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());

  1821 #else

  1822             if (e1 () (n, n) == value_type/*zero*/())

  1823                 singular ().raise ();

  1824 #endif

  1825             value_type t = e2 () (n) /= e1 () (n, n);

  1826             if (t != value_type/*zero*/()) {

  1827                 for (size_type m = n + 1; m < size; ++ m)

  1828                     e2 () (m) -= e1 () (m, n) * t;

  1829             }

  1830         }

  1831     }

  1832     // Packed (proxy) case

  1833     template<class E1, class E2>

  1834     BOOST_UBLAS_INLINE

  1835     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,

  1836                         lower_tag, column_major_tag, packed_proxy_tag) {

  1837         typedef typename E2::size_type size_type;

  1838         typedef typename E2::difference_type difference_type;

  1839         typedef typename E2::value_type value_type;

  1841         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());

  1842         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());

  1843         size_type size = e2 ().size ();

  1844         for (size_type n = 0; n < size; ++ n) {

  1845 #ifndef BOOST_UBLAS_SINGULAR_CHECK

  1846             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());

  1847 #else

  1848             if (e1 () (n, n) == value_type/*zero*/())

  1849                 singular ().raise ();

  1850 #endif

  1851             value_type t = e2 () (n) /= e1 () (n, n);

  1852             if (t != value_type/*zero*/()) {

  1853                 typename E1::const_iterator1 it1e1 (e1 ().find1 (1, n + 1, n));

  1854                 typename E1::const_iterator1 it1e1_end (e1 ().find1 (1, e1 ().size1 (), n));

  1855                 difference_type m (it1e1_end - it1e1);

  1856                 while (-- m >= 0)

  1857                     e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1;

  1858             }

  1859         }

  1860     }

  1861     // Sparse (proxy) case

  1862     template<class E1, class E2>

  1863     BOOST_UBLAS_INLINE

  1864     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,

  1865                         lower_tag, column_major_tag, unknown_storage_tag) {

  1866         typedef typename E2::size_type size_type;

  1867         typedef typename E2::difference_type difference_type;

  1868         typedef typename E2::value_type value_type;

  1870         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());

  1871         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());

  1872         size_type size = e2 ().size ();

  1873         for (size_type n = 0; n < size; ++ n) {

  1874 #ifndef BOOST_UBLAS_SINGULAR_CHECK

  1875             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());

  1876 #else

  1877             if (e1 () (n, n) == value_type/*zero*/())

  1878                 singular ().raise ();

  1879 #endif

  1880             value_type t = e2 () (n) /= e1 () (n, n);

  1881             if (t != value_type/*zero*/()) {

  1882                 typename E1::const_iterator1 it1e1 (e1 ().find1 (1, n + 1, n));

  1883                 typename E1::const_iterator1 it1e1_end (e1 ().find1 (1, e1 ().size1 (), n));

  1884                 while (it1e1 != it1e1_end)

  1885                     e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1;

  1886             }

  1887         }

  1888     }

  1889     // Redirectors :-)

  1890     template<class E1, class E2>

  1891     BOOST_UBLAS_INLINE

  1892     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,

  1893                         lower_tag, column_major_tag) {

  1894         typedef typename E1::storage_category storage_category;

  1895         inplace_solve (e1, e2,

  1896                        lower_tag (), column_major_tag (), storage_category ());

  1897     }

  1898     template<class E1, class E2>

  1899     BOOST_UBLAS_INLINE

  1900     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,

  1901                         lower_tag, row_major_tag) {

  1902         typedef typename E1::storage_category storage_category;

  1903         inplace_solve (e2, trans (e1),

  1904                        upper_tag (), row_major_tag (), storage_category ());

  1905     }

  1906     // Dispatcher

  1907     template<class E1, class E2>

  1908     BOOST_UBLAS_INLINE

  1909     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,

  1910                         lower_tag) {

  1911         typedef typename E1::orientation_category orientation_category;

  1912         inplace_solve (e1, e2,

  1913                        lower_tag (), orientation_category ());

  1914     }

  1915     template<class E1, class E2>

  1916     BOOST_UBLAS_INLINE

  1917     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,

  1918                         unit_lower_tag) {

  1919         typedef typename E1::orientation_category orientation_category;

  1920         inplace_solve (triangular_adaptor<const E1, unit_lower> (e1 ()), e2,

  1921                        unit_lower_tag (), orientation_category ());

  1922     }

  1924     // Dense (proxy) case

  1925     template<class E1, class E2>

  1926     BOOST_UBLAS_INLINE

  1927     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,

  1928                         upper_tag, column_major_tag, dense_proxy_tag) {

  1929         typedef typename E2::size_type size_type;

  1930         typedef typename E2::difference_type difference_type;

  1931         typedef typename E2::value_type value_type;

  1933         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());

  1934         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());

  1935         size_type size = e2 ().size ();

  1936         for (difference_type n = size - 1; n >= 0; -- n) {

  1937 #ifndef BOOST_UBLAS_SINGULAR_CHECK

  1938             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());

  1939 #else

  1940             if (e1 () (n, n) == value_type/*zero*/())

  1941                 singular ().raise ();

  1942 #endif

  1943             value_type t = e2 () (n) /= e1 () (n, n);

  1944             if (t != value_type/*zero*/()) {

  1945                 for (difference_type m = n - 1; m >= 0; -- m)

  1946                     e2 () (m) -= e1 () (m, n) * t;

  1947             }

  1948         }

  1949     }

  1950     // Packed (proxy) case

  1951     template<class E1, class E2>

  1952     BOOST_UBLAS_INLINE

  1953     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,

  1954                         upper_tag, column_major_tag, packed_proxy_tag) {

  1955         typedef typename E2::size_type size_type;

  1956         typedef typename E2::difference_type difference_type;

  1957         typedef typename E2::value_type value_type;

  1959         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());

  1960         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());

  1961         size_type size = e2 ().size ();

  1962         for (difference_type n = size - 1; n >= 0; -- n) {

  1963 #ifndef BOOST_UBLAS_SINGULAR_CHECK

  1964             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());

  1965 #else

  1966             if (e1 () (n, n) == value_type/*zero*/())

  1967                 singular ().raise ();

  1968 #endif

  1969             value_type t = e2 () (n) /= e1 () (n, n);

  1970             if (t != value_type/*zero*/()) {

  1971                 typename E1::const_reverse_iterator1 it1e1 (e1 ().find1 (1, n, n));

  1972                 typename E1::const_reverse_iterator1 it1e1_rend (e1 ().find1 (1, 0, n));

  1973                 difference_type m (it1e1_rend - it1e1);

  1974                 while (-- m >= 0)

  1975                     e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1;

  1976             }

  1977         }

  1978     }

  1979     // Sparse (proxy) case

  1980     template<class E1, class E2>

  1981     BOOST_UBLAS_INLINE

  1982     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,

  1983                         upper_tag, column_major_tag, unknown_storage_tag) {

  1984         typedef typename E2::size_type size_type;

  1985         typedef typename E2::difference_type difference_type;

  1986         typedef typename E2::value_type value_type;

  1988         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());

  1989         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());

  1990         size_type size = e2 ().size ();

  1991         for (difference_type n = size - 1; n >= 0; -- n) {

  1992 #ifndef BOOST_UBLAS_SINGULAR_CHECK

  1993             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());

  1994 #else

  1995             if (e1 () (n, n) == value_type/*zero*/())

  1996                 singular ().raise ();

  1997 #endif

  1998             value_type t = e2 () (n) /= e1 () (n, n);

  1999             if (t != value_type/*zero*/()) {

  2000                 typename E1::const_reverse_iterator1 it1e1 (e1 ().find1 (1, n, n));

  2001                 typename E1::const_reverse_iterator1 it1e1_rend (e1 ().find1 (1, 0, n));

  2002                 while (it1e1 != it1e1_rend)

  2003                     e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1;

  2004             }

  2005         }

  2006     }

  2007     // Redirectors :-)

  2008     template<class E1, class E2>

  2009     BOOST_UBLAS_INLINE

  2010     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,

  2011                         upper_tag, column_major_tag) {

  2012         typedef typename E1::storage_category storage_category;

  2013         inplace_solve (e1, e2,

  2014                        upper_tag (), column_major_tag (), storage_category ());

  2015     }

  2016     template<class E1, class E2>

  2017     BOOST_UBLAS_INLINE

  2018     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,

  2019                         upper_tag, row_major_tag) {

  2020         typedef typename E1::storage_category storage_category;

  2021         inplace_solve (e2, trans (e1),

  2022                        lower_tag (), row_major_tag (), storage_category ());

  2023     }

  2024     // Dispatcher

  2025     template<class E1, class E2>

  2026     BOOST_UBLAS_INLINE

  2027     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,

  2028                         upper_tag) {

  2029         typedef typename E1::orientation_category orientation_category;

  2030         inplace_solve (e1, e2,

  2031                        upper_tag (), orientation_category ());

  2032     }

  2033     template<class E1, class E2>

  2034     BOOST_UBLAS_INLINE

  2035     void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,

  2036                         unit_upper_tag) {

  2037         typedef typename E1::orientation_category orientation_category;

  2038         inplace_solve (triangular_adaptor<const E1, unit_upper> (e1 ()), e2,

  2039                        unit_upper_tag (), orientation_category ());

  2040     }

  2042     template<class E1, class E2, class C>

  2043     BOOST_UBLAS_INLINE

  2044     typename matrix_vector_solve_traits<E1, E2>::result_type

  2045     solve (const matrix_expression<E1> &e1,

  2046            const vector_expression<E2> &e2,

  2047            C) {

  2048         typename matrix_vector_solve_traits<E1, E2>::result_type r (e2);

  2049         inplace_solve (e1, r, C ());

  2050         return r;

  2051     }

  2053     // Dense (proxy) case

  2054     template<class E1, class E2>

  2055     BOOST_UBLAS_INLINE

  2056     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,

  2057                         lower_tag, row_major_tag, dense_proxy_tag) {

  2058         typedef typename E1::size_type size_type;

  2059         typedef typename E1::difference_type difference_type;

  2060         typedef typename E1::value_type value_type;

  2062         BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());

  2063         BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());

  2064         size_type size = e1 ().size ();

  2065         for (difference_type n = size - 1; n >= 0; -- n) {

  2066 #ifndef BOOST_UBLAS_SINGULAR_CHECK

  2067             BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());

  2068 #else

  2069             if (e2 () (n, n) == value_type/*zero*/())

  2070                 singular ().raise ();

  2071 #endif

  2072             value_type t = e1 () (n) /= e2 () (n, n);

  2073             if (t != value_type/*zero*/()) {

  2074                 for (difference_type m = n - 1; m >= 0; -- m)

  2075                     e1 () (m) -= t * e2 () (n, m);

  2076             }

  2077         }

  2078     }

  2079     // Packed (proxy) case

  2080     template<class E1, class E2>

  2081     BOOST_UBLAS_INLINE

  2082     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,

  2083                         lower_tag, row_major_tag, packed_proxy_tag) {

  2084         typedef typename E1::size_type size_type;

  2085         typedef typename E1::difference_type difference_type;

  2086         typedef typename E1::value_type value_type;

  2088         BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());

  2089         BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());

  2090         size_type size = e1 ().size ();

  2091         for (difference_type n = size - 1; n >= 0; -- n) {

  2092 #ifndef BOOST_UBLAS_SINGULAR_CHECK

  2093             BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());

  2094 #else

  2095             if (e2 () (n, n) == value_type/*zero*/())

  2096                 singular ().raise ();

  2097 #endif

  2098             value_type t = e1 () (n) /= e2 () (n, n);

  2099             if (t != value_type/*zero*/()) {

  2100                 typename E2::const_reverse_iterator2 it2e2 (e2 ().find2 (1, n, n));

  2101                 typename E2::const_reverse_iterator2 it2e2_rend (e2 ().find2 (1, n, 0));

  2102                 difference_type m (it2e2_rend - it2e2);

  2103                 while (-- m >= 0)

  2104                     e1 () (it2e2.index2 ()) -= *it2e2 * t, ++ it2e2;

  2105             }

  2106         }

  2107     }

  2108     // Sparse (proxy) case

  2109     template<class E1, class E2>

  2110     BOOST_UBLAS_INLINE

  2111     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,

  2112                         lower_tag, row_major_tag, unknown_storage_tag) {

  2113         typedef typename E1::size_type size_type;

  2114         typedef typename E1::difference_type difference_type;

  2115         typedef typename E1::value_type value_type;

  2117         BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());

  2118         BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());

  2119         size_type size = e1 ().size ();

  2120         for (difference_type n = size - 1; n >= 0; -- n) {

  2121 #ifndef BOOST_UBLAS_SINGULAR_CHECK

  2122             BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());

  2123 #else

  2124             if (e2 () (n, n) == value_type/*zero*/())

  2125                 singular ().raise ();

  2126 #endif

  2127             value_type t = e1 () (n) /= e2 () (n, n);

  2128             if (t != value_type/*zero*/()) {

  2129                 typename E2::const_reverse_iterator2 it2e2 (e2 ().find2 (1, n, n));

  2130                 typename E2::const_reverse_iterator2 it2e2_rend (e2 ().find2 (1, n, 0));

  2131                 while (it2e2 != it2e2_rend)

  2132                     e1 () (it2e2.index2 ()) -= *it2e2 * t, ++ it2e2;

  2133             }

  2134         }

  2135     }

  2136     // Redirectors :-)

  2137     template<class E1, class E2>

  2138     BOOST_UBLAS_INLINE

  2139     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,

  2140                         lower_tag, row_major_tag) {

  2141         typedef typename E1::storage_category storage_category;

  2142         inplace_solve (e1, e2,

  2143                        lower_tag (), row_major_tag (), storage_category ());

  2144     }

  2145     template<class E1, class E2>

  2146     BOOST_UBLAS_INLINE

  2147     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,

  2148                         lower_tag, column_major_tag) {

  2149         typedef typename E1::storage_category storage_category;

  2150         inplace_solve (trans (e2), e1,

  2151                        upper_tag (), row_major_tag (), storage_category ());

  2152     }

  2153     // Dispatcher

  2154     template<class E1, class E2>

  2155     BOOST_UBLAS_INLINE

  2156     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,

  2157                         lower_tag) {

  2158         typedef typename E2::orientation_category orientation_category;

  2159         inplace_solve (e1, e2,

  2160                        lower_tag (), orientation_category ());

  2161     }

  2162     template<class E1, class E2>

  2163     BOOST_UBLAS_INLINE

  2164     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,

  2165                         unit_lower_tag) {

  2166         typedef typename E2::orientation_category orientation_category;

  2167         inplace_solve (e1, triangular_adaptor<const E2, unit_lower> (e2 ()),

  2168                        unit_lower_tag (), orientation_category ());

  2169     }

  2171     // Dense (proxy) case

  2172     template<class E1, class E2>

  2173     BOOST_UBLAS_INLINE

  2174     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,

  2175                         upper_tag, row_major_tag, dense_proxy_tag) {

  2176         typedef typename E1::size_type size_type;

  2177         typedef typename E1::difference_type difference_type;

  2178         typedef typename E1::value_type value_type;

  2180         BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());

  2181         BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());

  2182         size_type size = e1 ().size ();

  2183         for (size_type n = 0; n < size; ++ n) {

  2184 #ifndef BOOST_UBLAS_SINGULAR_CHECK

  2185             BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());

  2186 #else

  2187             if (e2 () (n, n) == value_type/*zero*/())

  2188                 singular ().raise ();

  2189 #endif

  2190             value_type t = e1 () (n) /= e2 () (n, n);

  2191             if (t != value_type/*zero*/()) {

  2192                 for (size_type m = n + 1; m < size; ++ m)

  2193                     e1 () (m) -= t * e2 () (n, m);

  2194             }

  2195         }

  2196     }

  2197     // Packed (proxy) case

  2198     template<class E1, class E2>

  2199     BOOST_UBLAS_INLINE

  2200     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,

  2201                         upper_tag, row_major_tag, packed_proxy_tag) {

  2202         typedef typename E1::size_type size_type;

  2203         typedef typename E1::difference_type difference_type;

  2204         typedef typename E1::value_type value_type;

  2206         BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());

  2207         BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());

  2208         size_type size = e1 ().size ();

  2209         for (size_type n = 0; n < size; ++ n) {

  2210 #ifndef BOOST_UBLAS_SINGULAR_CHECK

  2211             BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());

  2212 #else

  2213             if (e2 () (n, n) == value_type/*zero*/())

  2214                 singular ().raise ();

  2215 #endif

  2216             value_type t = e1 () (n) /= e2 () (n, n);

  2217             if (t != value_type/*zero*/()) {

  2218                 typename E2::const_iterator2 it2e2 (e2 ().find2 (1, n, n + 1));

  2219                 typename E2::const_iterator2 it2e2_end (e2 ().find2 (1, n, e2 ().size2 ()));

  2220                 difference_type m (it2e2_end - it2e2);

  2221                 while (-- m >= 0)

  2222                     e1 () (it2e2.index2 ()) -= *it2e2 * t, ++ it2e2;

  2223             }

  2224         }

  2225     }

  2226     // Sparse (proxy) case

  2227     template<class E1, class E2>

  2228     BOOST_UBLAS_INLINE

  2229     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,

  2230                         upper_tag, row_major_tag, unknown_storage_tag) {

  2231         typedef typename E1::size_type size_type;

  2232         typedef typename E1::difference_type difference_type;

  2233         typedef typename E1::value_type value_type;

  2235         BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());

  2236         BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());

  2237         size_type size = e1 ().size ();

  2238         for (size_type n = 0; n < size; ++ n) {

  2239 #ifndef BOOST_UBLAS_SINGULAR_CHECK

  2240             BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());

  2241 #else

  2242             if (e2 () (n, n) == value_type/*zero*/())

  2243                 singular ().raise ();

  2244 #endif

  2245             value_type t = e1 () (n) /= e2 () (n, n);

  2246             if (t != value_type/*zero*/()) {

  2247                 typename E2::const_iterator2 it2e2 (e2 ().find2 (1, n, n + 1));

  2248                 typename E2::const_iterator2 it2e2_end (e2 ().find2 (1, n, e2 ().size2 ()));

  2249                 while (it2e2 != it2e2_end)

  2250                     e1 () (it2e2.index2 ()) -= *it2e2 * t, ++ it2e2;

  2251             }

  2252         }

  2253     }

  2254     // Redirectors :-)

  2255     template<class E1, class E2>

  2256     BOOST_UBLAS_INLINE

  2257     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,

  2258                         upper_tag, row_major_tag) {

  2259         typedef typename E1::storage_category storage_category;

  2260         inplace_solve (e1, e2,

  2261                        upper_tag (), row_major_tag (), storage_category ());

  2262     }

  2263     template<class E1, class E2>

  2264     BOOST_UBLAS_INLINE

  2265     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,

  2266                         upper_tag, column_major_tag) {

  2267         typedef typename E1::storage_category storage_category;

  2268         inplace_solve (trans (e2), e1,

  2269                        lower_tag (), row_major_tag (), storage_category ());

  2270     }

  2271     // Dispatcher

  2272     template<class E1, class E2>

  2273     BOOST_UBLAS_INLINE

  2274     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,

  2275                         upper_tag) {

  2276         typedef typename E2::orientation_category orientation_category;

  2277         inplace_solve (e1, e2,

  2278                        upper_tag (), orientation_category ());

  2279     }

  2280     template<class E1, class E2>

  2281     BOOST_UBLAS_INLINE

  2282     void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,

  2283                         unit_upper_tag) {

  2284         typedef typename E2::orientation_category orientation_category;

  2285         inplace_solve (e1, triangular_adaptor<const E2, unit_upper> (e2 ()),

  2286                        unit_upper_tag (), orientation_category ());

  2287     }

  2289     template<class E1, class E2, class C>

  2290     BOOST_UBLAS_INLINE

  2291     typename matrix_vector_solve_traits<E1, E2>::result_type

  2292     solve (const vector_expression<E1> &e1,

  2293            const matrix_expression<E2> &e2,

  2294            C) {

  2295         typename matrix_vector_solve_traits<E1, E2>::result_type r (e1);

  2296         inplace_solve (r, e2, C ());

  2297         return r;

  2298     }

  2300     template<class E1, class E2>

  2301     struct matrix_matrix_solve_traits {

  2302         typedef typename promote_traits<typename E1::value_type, typename E2::value_type>::promote_type promote_type;

  2303         typedef matrix<promote_type> result_type;

  2304     };

  2306     // Operations:

  2307     //  k * n * (n - 1) / 2 + k * n = k * n * (n + 1) / 2 multiplications,

  2308     //  k * n * (n - 1) / 2 additions

  2310     // Dense (proxy) case

  2311     template<class E1, class E2>

  2312     BOOST_UBLAS_INLINE

  2313     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,

  2314                         lower_tag, dense_proxy_tag) {

  2315         typedef typename E2::size_type size_type;

  2316         typedef typename E2::difference_type difference_type;

  2317         typedef typename E2::value_type value_type;

  2319         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());

  2320         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());

  2321         size_type size1 = e2 ().size1 ();

  2322         size_type size2 = e2 ().size2 ();

  2323         for (size_type n = 0; n < size1; ++ n) {

  2324 #ifndef BOOST_UBLAS_SINGULAR_CHECK

  2325             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());

  2326 #else

  2327             if (e1 () (n, n) == value_type/*zero*/())

  2328                 singular ().raise ();

  2329 #endif

  2330             for (size_type l = 0; l < size2; ++ l) {

  2331                 value_type t = e2 () (n, l) /= e1 () (n, n);

  2332                 if (t != value_type/*zero*/()) {

  2333                     for (size_type m = n + 1; m < size1; ++ m)

  2334                         e2 () (m, l) -= e1 () (m, n) * t;

  2335                 }

  2336             }

  2337         }

  2338     }

  2339     // Packed (proxy) case

  2340     template<class E1, class E2>

  2341     BOOST_UBLAS_INLINE

  2342     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,

  2343                         lower_tag, packed_proxy_tag) {

  2344         typedef typename E2::size_type size_type;

  2345         typedef typename E2::difference_type difference_type;

  2346         typedef typename E2::value_type value_type;

  2348         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());

  2349         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());

  2350         size_type size1 = e2 ().size1 ();

  2351         size_type size2 = e2 ().size2 ();

  2352         for (size_type n = 0; n < size1; ++ n) {

  2353 #ifndef BOOST_UBLAS_SINGULAR_CHECK

  2354             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());

  2355 #else

  2356             if (e1 () (n, n) == value_type/*zero*/())

  2357                 singular ().raise ();

  2358 #endif

  2359             for (size_type l = 0; l < size2; ++ l) {

  2360                 value_type t = e2 () (n, l) /= e1 () (n, n);

  2361                 if (t != value_type/*zero*/()) {

  2362                     typename E1::const_iterator1 it1e1 (e1 ().find1 (1, n + 1, n));

  2363                     typename E1::const_iterator1 it1e1_end (e1 ().find1 (1, e1 ().size1 (), n));

  2364                     difference_type m (it1e1_end - it1e1);

  2365                     while (-- m >= 0)

  2366                         e2 () (it1e1.index1 (), l) -= *it1e1 * t, ++ it1e1;

  2367                 }

  2368             }

  2369         }

  2370     }

  2371     // Sparse (proxy) case

  2372     template<class E1, class E2>

  2373     BOOST_UBLAS_INLINE

  2374     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,

  2375                         lower_tag, unknown_storage_tag) {

  2376         typedef typename E2::size_type size_type;

  2377         typedef typename E2::difference_type difference_type;

  2378         typedef typename E2::value_type value_type;

  2380         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());

  2381         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());

  2382         size_type size1 = e2 ().size1 ();

  2383         size_type size2 = e2 ().size2 ();

  2384         for (size_type n = 0; n < size1; ++ n) {

  2385 #ifndef BOOST_UBLAS_SINGULAR_CHECK

  2386             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());

  2387 #else

  2388             if (e1 () (n, n) == value_type/*zero*/())

  2389                 singular ().raise ();

  2390 #endif

  2391             for (size_type l = 0; l < size2; ++ l) {

  2392                 value_type t = e2 () (n, l) /= e1 () (n, n);

  2393                 if (t != value_type/*zero*/()) {

  2394                     typename E1::const_iterator1 it1e1 (e1 ().find1 (1, n + 1, n));

  2395                     typename E1::const_iterator1 it1e1_end (e1 ().find1 (1, e1 ().size1 (), n));

  2396                     while (it1e1 != it1e1_end)

  2397                         e2 () (it1e1.index1 (), l) -= *it1e1 * t, ++ it1e1;

  2398                 }

  2399             }

  2400         }

  2401     }

  2402     // Dispatcher

  2403     template<class E1, class E2>

  2404     BOOST_UBLAS_INLINE

  2405     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,

  2406                         lower_tag) {

  2407         typedef typename E1::storage_category dispatch_category;

  2408         inplace_solve (e1, e2,

  2409                        lower_tag (), dispatch_category ());

  2410     }

  2411     template<class E1, class E2>

  2412     BOOST_UBLAS_INLINE

  2413     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,

  2414                         unit_lower_tag) {

  2415         typedef typename E1::storage_category dispatch_category;

  2416         inplace_solve (triangular_adaptor<const E1, unit_lower> (e1 ()), e2,

  2417                        unit_lower_tag (), dispatch_category ());

  2418     }

  2420     // Dense (proxy) case

  2421     template<class E1, class E2>

  2422     BOOST_UBLAS_INLINE

  2423     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,

  2424                         upper_tag, dense_proxy_tag) {

  2425         typedef typename E2::size_type size_type;

  2426         typedef typename E2::difference_type difference_type;

  2427         typedef typename E2::value_type value_type;

  2429         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());

  2430         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());

  2431         size_type size1 = e2 ().size1 ();

  2432         size_type size2 = e2 ().size2 ();

  2433         for (difference_type n = size1 - 1; n >= 0; -- n) {

  2434 #ifndef BOOST_UBLAS_SINGULAR_CHECK

  2435             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());

  2436 #else

  2437             if (e1 () (n, n) == value_type/*zero*/())

  2438                 singular ().raise ();

  2439 #endif

  2440             for (difference_type l = size2 - 1; l >= 0; -- l) {

  2441                 value_type t = e2 () (n, l) /= e1 () (n, n);

  2442                 if (t != value_type/*zero*/()) {

  2443                     for (difference_type m = n - 1; m >= 0; -- m)

  2444                         e2 () (m, l) -= e1 () (m, n) * t;

  2445                 }

  2446             }

  2447         }

  2448     }

  2449     // Packed (proxy) case

  2450     template<class E1, class E2>

  2451     BOOST_UBLAS_INLINE

  2452     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,

  2453                         upper_tag, packed_proxy_tag) {

  2454         typedef typename E2::size_type size_type;

  2455         typedef typename E2::difference_type difference_type;

  2456         typedef typename E2::value_type value_type;

  2458         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());

  2459         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());

  2460         size_type size1 = e2 ().size1 ();

  2461         size_type size2 = e2 ().size2 ();

  2462         for (difference_type n = size1 - 1; n >= 0; -- n) {

  2463 #ifndef BOOST_UBLAS_SINGULAR_CHECK

  2464             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());

  2465 #else

  2466             if (e1 () (n, n) == value_type/*zero*/())

  2467                 singular ().raise ();

  2468 #endif

  2469             for (difference_type l = size2 - 1; l >= 0; -- l) {

  2470                 value_type t = e2 () (n, l) /= e1 () (n, n);

  2471                 if (t != value_type/*zero*/()) {

  2472                     typename E1::const_reverse_iterator1 it1e1 (e1 ().find1 (1, n, n));

  2473                     typename E1::const_reverse_iterator1 it1e1_rend (e1 ().find1 (1, 0, n));

  2474                     difference_type m (it1e1_rend - it1e1);

  2475                     while (-- m >= 0)

  2476                         e2 () (it1e1.index1 (), l) -= *it1e1 * t, ++ it1e1;

  2477                 }

  2478             }

  2479         }

  2480     }

  2481     // Sparse (proxy) case

  2482     template<class E1, class E2>

  2483     BOOST_UBLAS_INLINE

  2484     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,

  2485                         upper_tag, unknown_storage_tag) {

  2486         typedef typename E2::size_type size_type;

  2487         typedef typename E2::difference_type difference_type;

  2488         typedef typename E2::value_type value_type;

  2490         BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());

  2491         BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());

  2492         size_type size1 = e2 ().size1 ();

  2493         size_type size2 = e2 ().size2 ();

  2494         for (difference_type n = size1 - 1; n >= 0; -- n) {

  2495 #ifndef BOOST_UBLAS_SINGULAR_CHECK

  2496             BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());

  2497 #else

  2498             if (e1 () (n, n) == value_type/*zero*/())

  2499                 singular ().raise ();

  2500 #endif

  2501             for (difference_type l = size2 - 1; l >= 0; -- l) {

  2502                 value_type t = e2 () (n, l) /= e1 () (n, n);

  2503                 if (t != value_type/*zero*/()) {

  2504                     typename E1::const_reverse_iterator1 it1e1 (e1 ().find1 (1, n, n));

  2505                     typename E1::const_reverse_iterator1 it1e1_rend (e1 ().find1 (1, 0, n));

  2506                     while (it1e1 != it1e1_rend)

  2507                         e2 () (it1e1.index1 (), l) -= *it1e1 * t, ++ it1e1;

  2508                 }

  2509             }

  2510         }

  2511     }

  2512     // Dispatcher

  2513     template<class E1, class E2>

  2514     BOOST_UBLAS_INLINE

  2515     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,

  2516                         upper_tag) {

  2517         typedef typename E1::storage_category dispatch_category;

  2518         inplace_solve (e1, e2,

  2519                        upper_tag (), dispatch_category ());

  2520     }

  2521     template<class E1, class E2>

  2522     BOOST_UBLAS_INLINE

  2523     void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,

  2524                         unit_upper_tag) {

  2525         typedef typename E1::storage_category dispatch_category;

  2526         inplace_solve (triangular_adaptor<const E1, unit_upper> (e1 ()), e2,

  2527                        unit_upper_tag (), dispatch_category ());

  2528     }

  2530     template<class E1, class E2, class C>

  2531     BOOST_UBLAS_INLINE

  2532     typename matrix_matrix_solve_traits<E1, E2>::result_type

  2533     solve (const matrix_expression<E1> &e1,

  2534            const matrix_expression<E2> &e2,

  2535            C) {

  2536         typename matrix_matrix_solve_traits<E1, E2>::result_type r (e2);

  2537         inplace_solve (e1, r, C ());

  2538         return r;

  2539     }

  2541 }}}

  2543 #endif

author	sl
	Tue, 10 Jun 2014 14:32:02 +0200
changeset 1	260cb5ec6c19
permissions	-rw-r--r--