1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/os/ossrv/ossrv_pub/boost_apis/boost/numeric/ublas/triangular.hpp Fri Jun 15 03:10:57 2012 +0200
1.3 @@ -0,0 +1,2543 @@
1.4 +//
1.5 +// Copyright (c) 2000-2002
1.6 +// Joerg Walter, Mathias Koch
1.7 +//
1.8 +// Permission to use, copy, modify, distribute and sell this software
1.9 +// and its documentation for any purpose is hereby granted without fee,
1.10 +// provided that the above copyright notice appear in all copies and
1.11 +// that both that copyright notice and this permission notice appear
1.12 +// in supporting documentation. The authors make no representations
1.13 +// about the suitability of this software for any purpose.
1.14 +// It is provided "as is" without express or implied warranty.
1.15 +//
1.16 +// The authors gratefully acknowledge the support of
1.17 +// GeNeSys mbH & Co. KG in producing this work.
1.18 +//
1.19 +
1.20 +#ifndef _BOOST_UBLAS_TRIANGULAR_
1.21 +#define _BOOST_UBLAS_TRIANGULAR_
1.22 +
1.23 +#include <boost/numeric/ublas/matrix.hpp>
1.24 +#include <boost/numeric/ublas/detail/temporary.hpp>
1.25 +#include <boost/type_traits/remove_const.hpp>
1.26 +
1.27 +// Iterators based on ideas of Jeremy Siek
1.28 +
1.29 +namespace boost { namespace numeric { namespace ublas {
1.30 +
1.31 + // Array based triangular matrix class
1.32 + template<class T, class TRI, class L, class A>
1.33 + class triangular_matrix:
1.34 + public matrix_container<triangular_matrix<T, TRI, L, A> > {
1.35 +
1.36 + typedef T *pointer;
1.37 + typedef TRI triangular_type;
1.38 + typedef L layout_type;
1.39 + typedef triangular_matrix<T, TRI, L, A> self_type;
1.40 + public:
1.41 +#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
1.42 + using matrix_container<self_type>::operator ();
1.43 +#endif
1.44 + typedef typename A::size_type size_type;
1.45 + typedef typename A::difference_type difference_type;
1.46 + typedef T value_type;
1.47 + typedef const T &const_reference;
1.48 + typedef T &reference;
1.49 + typedef A array_type;
1.50 +
1.51 + typedef const matrix_reference<const self_type> const_closure_type;
1.52 + typedef matrix_reference<self_type> closure_type;
1.53 + typedef vector<T, A> vector_temporary_type;
1.54 + typedef matrix<T, L, A> matrix_temporary_type; // general sub-matrix
1.55 + typedef packed_tag storage_category;
1.56 + typedef typename L::orientation_category orientation_category;
1.57 +
1.58 + // Construction and destruction
1.59 + BOOST_UBLAS_INLINE
1.60 + triangular_matrix ():
1.61 + matrix_container<self_type> (),
1.62 + size1_ (0), size2_ (0), data_ (0) {}
1.63 + BOOST_UBLAS_INLINE
1.64 + triangular_matrix (size_type size1, size_type size2):
1.65 + matrix_container<self_type> (),
1.66 + size1_ (size1), size2_ (size2), data_ (triangular_type::packed_size (layout_type (), size1, size2)) {
1.67 + }
1.68 + BOOST_UBLAS_INLINE
1.69 + triangular_matrix (size_type size1, size_type size2, const array_type &data):
1.70 + matrix_container<self_type> (),
1.71 + size1_ (size1), size2_ (size2), data_ (data) {}
1.72 + BOOST_UBLAS_INLINE
1.73 + triangular_matrix (const triangular_matrix &m):
1.74 + matrix_container<self_type> (),
1.75 + size1_ (m.size1_), size2_ (m.size2_), data_ (m.data_) {}
1.76 + template<class AE>
1.77 + BOOST_UBLAS_INLINE
1.78 + triangular_matrix (const matrix_expression<AE> &ae):
1.79 + matrix_container<self_type> (),
1.80 + size1_ (ae ().size1 ()), size2_ (ae ().size2 ()),
1.81 + data_ (triangular_type::packed_size (layout_type (), size1_, size2_)) {
1.82 + matrix_assign<scalar_assign> (*this, ae);
1.83 + }
1.84 +
1.85 + // Accessors
1.86 + BOOST_UBLAS_INLINE
1.87 + size_type size1 () const {
1.88 + return size1_;
1.89 + }
1.90 + BOOST_UBLAS_INLINE
1.91 + size_type size2 () const {
1.92 + return size2_;
1.93 + }
1.94 +
1.95 + // Storage accessors
1.96 + BOOST_UBLAS_INLINE
1.97 + const array_type &data () const {
1.98 + return data_;
1.99 + }
1.100 + BOOST_UBLAS_INLINE
1.101 + array_type &data () {
1.102 + return data_;
1.103 + }
1.104 +
1.105 + // Resizing
1.106 + BOOST_UBLAS_INLINE
1.107 + void resize (size_type size1, size_type size2, bool preserve = true) {
1.108 + if (preserve) {
1.109 + self_type temporary (size1, size2);
1.110 + detail::matrix_resize_preserve<layout_type> (*this, temporary);
1.111 + }
1.112 + else {
1.113 + data ().resize (triangular_type::packed_size (layout_type (), size1, size2));
1.114 + size1_ = size1;
1.115 + size2_ = size2;
1.116 + }
1.117 + }
1.118 + BOOST_UBLAS_INLINE
1.119 + void resize_packed_preserve (size_type size1, size_type size2) {
1.120 + size1_ = size1;
1.121 + size2_ = size2;
1.122 + data ().resize (triangular_type::packed_size (layout_type (), size1_, size2_), value_type ());
1.123 + }
1.124 +
1.125 + // Element access
1.126 + BOOST_UBLAS_INLINE
1.127 + const_reference operator () (size_type i, size_type j) const {
1.128 + BOOST_UBLAS_CHECK (i < size1_, bad_index ());
1.129 + BOOST_UBLAS_CHECK (j < size2_, bad_index ());
1.130 + if (triangular_type::other (i, j))
1.131 + return data () [triangular_type::element (layout_type (), i, size1_, j, size2_)];
1.132 + else if (triangular_type::one (i, j))
1.133 + return one_;
1.134 + else
1.135 + return zero_;
1.136 + }
1.137 + BOOST_UBLAS_INLINE
1.138 + reference at_element (size_type i, size_type j) {
1.139 + BOOST_UBLAS_CHECK (i < size1_, bad_index ());
1.140 + BOOST_UBLAS_CHECK (j < size2_, bad_index ());
1.141 + return data () [triangular_type::element (layout_type (), i, size1_, j, size2_)];
1.142 + }
1.143 + BOOST_UBLAS_INLINE
1.144 + reference operator () (size_type i, size_type j) {
1.145 + BOOST_UBLAS_CHECK (i < size1_, bad_index ());
1.146 + BOOST_UBLAS_CHECK (j < size2_, bad_index ());
1.147 + if (triangular_type::other (i, j))
1.148 + return data () [triangular_type::element (layout_type (), i, size1_, j, size2_)];
1.149 + else {
1.150 + bad_index ().raise ();
1.151 + // arbitary return value
1.152 + return const_cast<reference>(zero_);
1.153 + }
1.154 + }
1.155 +
1.156 + // Element assignment
1.157 + BOOST_UBLAS_INLINE
1.158 + reference insert_element (size_type i, size_type j, const_reference t) {
1.159 + return (operator () (i, j) = t);
1.160 + }
1.161 + BOOST_UBLAS_INLINE
1.162 + void erase_element (size_type i, size_type j) {
1.163 + operator () (i, j) = value_type/*zero*/();
1.164 + }
1.165 +
1.166 + // Zeroing
1.167 + BOOST_UBLAS_INLINE
1.168 + void clear () {
1.169 + // data ().clear ();
1.170 + std::fill (data ().begin (), data ().end (), value_type/*zero*/());
1.171 + }
1.172 +
1.173 + // Assignment
1.174 + BOOST_UBLAS_INLINE
1.175 + triangular_matrix &operator = (const triangular_matrix &m) {
1.176 + size1_ = m.size1_;
1.177 + size2_ = m.size2_;
1.178 + data () = m.data ();
1.179 + return *this;
1.180 + }
1.181 + BOOST_UBLAS_INLINE
1.182 + triangular_matrix &assign_temporary (triangular_matrix &m) {
1.183 + swap (m);
1.184 + return *this;
1.185 + }
1.186 + template<class AE>
1.187 + BOOST_UBLAS_INLINE
1.188 + triangular_matrix &operator = (const matrix_expression<AE> &ae) {
1.189 + self_type temporary (ae);
1.190 + return assign_temporary (temporary);
1.191 + }
1.192 + template<class AE>
1.193 + BOOST_UBLAS_INLINE
1.194 + triangular_matrix &assign (const matrix_expression<AE> &ae) {
1.195 + matrix_assign<scalar_assign> (*this, ae);
1.196 + return *this;
1.197 + }
1.198 + template<class AE>
1.199 + BOOST_UBLAS_INLINE
1.200 + triangular_matrix& operator += (const matrix_expression<AE> &ae) {
1.201 + self_type temporary (*this + ae);
1.202 + return assign_temporary (temporary);
1.203 + }
1.204 + template<class AE>
1.205 + BOOST_UBLAS_INLINE
1.206 + triangular_matrix &plus_assign (const matrix_expression<AE> &ae) {
1.207 + matrix_assign<scalar_plus_assign> (*this, ae);
1.208 + return *this;
1.209 + }
1.210 + template<class AE>
1.211 + BOOST_UBLAS_INLINE
1.212 + triangular_matrix& operator -= (const matrix_expression<AE> &ae) {
1.213 + self_type temporary (*this - ae);
1.214 + return assign_temporary (temporary);
1.215 + }
1.216 + template<class AE>
1.217 + BOOST_UBLAS_INLINE
1.218 + triangular_matrix &minus_assign (const matrix_expression<AE> &ae) {
1.219 + matrix_assign<scalar_minus_assign> (*this, ae);
1.220 + return *this;
1.221 + }
1.222 + template<class AT>
1.223 + BOOST_UBLAS_INLINE
1.224 + triangular_matrix& operator *= (const AT &at) {
1.225 + matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
1.226 + return *this;
1.227 + }
1.228 + template<class AT>
1.229 + BOOST_UBLAS_INLINE
1.230 + triangular_matrix& operator /= (const AT &at) {
1.231 + matrix_assign_scalar<scalar_divides_assign> (*this, at);
1.232 + return *this;
1.233 + }
1.234 +
1.235 + // Swapping
1.236 + BOOST_UBLAS_INLINE
1.237 + void swap (triangular_matrix &m) {
1.238 + if (this != &m) {
1.239 + // BOOST_UBLAS_CHECK (size2_ == m.size2_, bad_size ());
1.240 + std::swap (size1_, m.size1_);
1.241 + std::swap (size2_, m.size2_);
1.242 + data ().swap (m.data ());
1.243 + }
1.244 + }
1.245 + BOOST_UBLAS_INLINE
1.246 + friend void swap (triangular_matrix &m1, triangular_matrix &m2) {
1.247 + m1.swap (m2);
1.248 + }
1.249 +
1.250 + // Iterator types
1.251 +#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1.252 + typedef indexed_iterator1<self_type, packed_random_access_iterator_tag> iterator1;
1.253 + typedef indexed_iterator2<self_type, packed_random_access_iterator_tag> iterator2;
1.254 + typedef indexed_const_iterator1<self_type, packed_random_access_iterator_tag> const_iterator1;
1.255 + typedef indexed_const_iterator2<self_type, packed_random_access_iterator_tag> const_iterator2;
1.256 +#else
1.257 + class const_iterator1;
1.258 + class iterator1;
1.259 + class const_iterator2;
1.260 + class iterator2;
1.261 +#endif
1.262 + typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
1.263 + typedef reverse_iterator_base1<iterator1> reverse_iterator1;
1.264 + typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
1.265 + typedef reverse_iterator_base2<iterator2> reverse_iterator2;
1.266 +
1.267 + // Element lookup
1.268 + BOOST_UBLAS_INLINE
1.269 + const_iterator1 find1 (int rank, size_type i, size_type j) const {
1.270 + if (rank == 1)
1.271 + i = triangular_type::restrict1 (i, j);
1.272 + return const_iterator1 (*this, i, j);
1.273 + }
1.274 + BOOST_UBLAS_INLINE
1.275 + iterator1 find1 (int rank, size_type i, size_type j) {
1.276 + if (rank == 1)
1.277 + i = triangular_type::mutable_restrict1 (i, j);
1.278 + return iterator1 (*this, i, j);
1.279 + }
1.280 + BOOST_UBLAS_INLINE
1.281 + const_iterator2 find2 (int rank, size_type i, size_type j) const {
1.282 + if (rank == 1)
1.283 + j = triangular_type::restrict2 (i, j);
1.284 + return const_iterator2 (*this, i, j);
1.285 + }
1.286 + BOOST_UBLAS_INLINE
1.287 + iterator2 find2 (int rank, size_type i, size_type j) {
1.288 + if (rank == 1)
1.289 + j = triangular_type::mutable_restrict2 (i, j);
1.290 + return iterator2 (*this, i, j);
1.291 + }
1.292 +
1.293 + // Iterators simply are indices.
1.294 +
1.295 +#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1.296 + class const_iterator1:
1.297 + public container_const_reference<triangular_matrix>,
1.298 + public random_access_iterator_base<packed_random_access_iterator_tag,
1.299 + const_iterator1, value_type> {
1.300 + public:
1.301 + typedef typename triangular_matrix::value_type value_type;
1.302 + typedef typename triangular_matrix::difference_type difference_type;
1.303 + typedef typename triangular_matrix::const_reference reference;
1.304 + typedef const typename triangular_matrix::pointer pointer;
1.305 +
1.306 + typedef const_iterator2 dual_iterator_type;
1.307 + typedef const_reverse_iterator2 dual_reverse_iterator_type;
1.308 +
1.309 + // Construction and destruction
1.310 + BOOST_UBLAS_INLINE
1.311 + const_iterator1 ():
1.312 + container_const_reference<self_type> (), it1_ (), it2_ () {}
1.313 + BOOST_UBLAS_INLINE
1.314 + const_iterator1 (const self_type &m, size_type it1, size_type it2):
1.315 + container_const_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
1.316 + BOOST_UBLAS_INLINE
1.317 + const_iterator1 (const iterator1 &it):
1.318 + container_const_reference<self_type> (it ()), it1_ (it.it1_), it2_ (it.it2_) {}
1.319 +
1.320 + // Arithmetic
1.321 + BOOST_UBLAS_INLINE
1.322 + const_iterator1 &operator ++ () {
1.323 + ++ it1_;
1.324 + return *this;
1.325 + }
1.326 + BOOST_UBLAS_INLINE
1.327 + const_iterator1 &operator -- () {
1.328 + -- it1_;
1.329 + return *this;
1.330 + }
1.331 + BOOST_UBLAS_INLINE
1.332 + const_iterator1 &operator += (difference_type n) {
1.333 + it1_ += n;
1.334 + return *this;
1.335 + }
1.336 + BOOST_UBLAS_INLINE
1.337 + const_iterator1 &operator -= (difference_type n) {
1.338 + it1_ -= n;
1.339 + return *this;
1.340 + }
1.341 + BOOST_UBLAS_INLINE
1.342 + difference_type operator - (const const_iterator1 &it) const {
1.343 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.344 + BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
1.345 + return it1_ - it.it1_;
1.346 + }
1.347 +
1.348 + // Dereference
1.349 + BOOST_UBLAS_INLINE
1.350 + const_reference operator * () const {
1.351 + return (*this) () (it1_, it2_);
1.352 + }
1.353 + BOOST_UBLAS_INLINE
1.354 + const_reference operator [] (difference_type n) const {
1.355 + return *(*this + n);
1.356 + }
1.357 +
1.358 +#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1.359 + BOOST_UBLAS_INLINE
1.360 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.361 + typename self_type::
1.362 +#endif
1.363 + const_iterator2 begin () const {
1.364 + return (*this) ().find2 (1, it1_, 0);
1.365 + }
1.366 + BOOST_UBLAS_INLINE
1.367 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.368 + typename self_type::
1.369 +#endif
1.370 + const_iterator2 end () const {
1.371 + return (*this) ().find2 (1, it1_, (*this) ().size2 ());
1.372 + }
1.373 + BOOST_UBLAS_INLINE
1.374 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.375 + typename self_type::
1.376 +#endif
1.377 + const_reverse_iterator2 rbegin () const {
1.378 + return const_reverse_iterator2 (end ());
1.379 + }
1.380 + BOOST_UBLAS_INLINE
1.381 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.382 + typename self_type::
1.383 +#endif
1.384 + const_reverse_iterator2 rend () const {
1.385 + return const_reverse_iterator2 (begin ());
1.386 + }
1.387 +#endif
1.388 +
1.389 + // Indices
1.390 + BOOST_UBLAS_INLINE
1.391 + size_type index1 () const {
1.392 + return it1_;
1.393 + }
1.394 + BOOST_UBLAS_INLINE
1.395 + size_type index2 () const {
1.396 + return it2_;
1.397 + }
1.398 +
1.399 + // Assignment
1.400 + BOOST_UBLAS_INLINE
1.401 + const_iterator1 &operator = (const const_iterator1 &it) {
1.402 + container_const_reference<self_type>::assign (&it ());
1.403 + it1_ = it.it1_;
1.404 + it2_ = it.it2_;
1.405 + return *this;
1.406 + }
1.407 +
1.408 + // Comparison
1.409 + BOOST_UBLAS_INLINE
1.410 + bool operator == (const const_iterator1 &it) const {
1.411 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.412 + BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
1.413 + return it1_ == it.it1_;
1.414 + }
1.415 + BOOST_UBLAS_INLINE
1.416 + bool operator < (const const_iterator1 &it) const {
1.417 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.418 + BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
1.419 + return it1_ < it.it1_;
1.420 + }
1.421 +
1.422 + private:
1.423 + size_type it1_;
1.424 + size_type it2_;
1.425 + };
1.426 +#endif
1.427 +
1.428 + BOOST_UBLAS_INLINE
1.429 + const_iterator1 begin1 () const {
1.430 + return find1 (0, 0, 0);
1.431 + }
1.432 + BOOST_UBLAS_INLINE
1.433 + const_iterator1 end1 () const {
1.434 + return find1 (0, size1_, 0);
1.435 + }
1.436 +
1.437 +#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1.438 + class iterator1:
1.439 + public container_reference<triangular_matrix>,
1.440 + public random_access_iterator_base<packed_random_access_iterator_tag,
1.441 + iterator1, value_type> {
1.442 + public:
1.443 + typedef typename triangular_matrix::value_type value_type;
1.444 + typedef typename triangular_matrix::difference_type difference_type;
1.445 + typedef typename triangular_matrix::reference reference;
1.446 + typedef typename triangular_matrix::pointer pointer;
1.447 +
1.448 + typedef iterator2 dual_iterator_type;
1.449 + typedef reverse_iterator2 dual_reverse_iterator_type;
1.450 +
1.451 + // Construction and destruction
1.452 + BOOST_UBLAS_INLINE
1.453 + iterator1 ():
1.454 + container_reference<self_type> (), it1_ (), it2_ () {}
1.455 + BOOST_UBLAS_INLINE
1.456 + iterator1 (self_type &m, size_type it1, size_type it2):
1.457 + container_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
1.458 +
1.459 + // Arithmetic
1.460 + BOOST_UBLAS_INLINE
1.461 + iterator1 &operator ++ () {
1.462 + ++ it1_;
1.463 + return *this;
1.464 + }
1.465 + BOOST_UBLAS_INLINE
1.466 + iterator1 &operator -- () {
1.467 + -- it1_;
1.468 + return *this;
1.469 + }
1.470 + BOOST_UBLAS_INLINE
1.471 + iterator1 &operator += (difference_type n) {
1.472 + it1_ += n;
1.473 + return *this;
1.474 + }
1.475 + BOOST_UBLAS_INLINE
1.476 + iterator1 &operator -= (difference_type n) {
1.477 + it1_ -= n;
1.478 + return *this;
1.479 + }
1.480 + BOOST_UBLAS_INLINE
1.481 + difference_type operator - (const iterator1 &it) const {
1.482 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.483 + BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
1.484 + return it1_ - it.it1_;
1.485 + }
1.486 +
1.487 + // Dereference
1.488 + BOOST_UBLAS_INLINE
1.489 + reference operator * () const {
1.490 + return (*this) () (it1_, it2_);
1.491 + }
1.492 + BOOST_UBLAS_INLINE
1.493 + reference operator [] (difference_type n) const {
1.494 + return *(*this + n);
1.495 + }
1.496 +
1.497 +#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1.498 + BOOST_UBLAS_INLINE
1.499 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.500 + typename self_type::
1.501 +#endif
1.502 + iterator2 begin () const {
1.503 + return (*this) ().find2 (1, it1_, 0);
1.504 + }
1.505 + BOOST_UBLAS_INLINE
1.506 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.507 + typename self_type::
1.508 +#endif
1.509 + iterator2 end () const {
1.510 + return (*this) ().find2 (1, it1_, (*this) ().size2 ());
1.511 + }
1.512 + BOOST_UBLAS_INLINE
1.513 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.514 + typename self_type::
1.515 +#endif
1.516 + reverse_iterator2 rbegin () const {
1.517 + return reverse_iterator2 (end ());
1.518 + }
1.519 + BOOST_UBLAS_INLINE
1.520 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.521 + typename self_type::
1.522 +#endif
1.523 + reverse_iterator2 rend () const {
1.524 + return reverse_iterator2 (begin ());
1.525 + }
1.526 +#endif
1.527 +
1.528 + // Indices
1.529 + BOOST_UBLAS_INLINE
1.530 + size_type index1 () const {
1.531 + return it1_;
1.532 + }
1.533 + BOOST_UBLAS_INLINE
1.534 + size_type index2 () const {
1.535 + return it2_;
1.536 + }
1.537 +
1.538 + // Assignment
1.539 + BOOST_UBLAS_INLINE
1.540 + iterator1 &operator = (const iterator1 &it) {
1.541 + container_reference<self_type>::assign (&it ());
1.542 + it1_ = it.it1_;
1.543 + it2_ = it.it2_;
1.544 + return *this;
1.545 + }
1.546 +
1.547 + // Comparison
1.548 + BOOST_UBLAS_INLINE
1.549 + bool operator == (const iterator1 &it) const {
1.550 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.551 + BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
1.552 + return it1_ == it.it1_;
1.553 + }
1.554 + BOOST_UBLAS_INLINE
1.555 + bool operator < (const iterator1 &it) const {
1.556 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.557 + BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ());
1.558 + return it1_ < it.it1_;
1.559 + }
1.560 +
1.561 + private:
1.562 + size_type it1_;
1.563 + size_type it2_;
1.564 +
1.565 + friend class const_iterator1;
1.566 + };
1.567 +#endif
1.568 +
1.569 + BOOST_UBLAS_INLINE
1.570 + iterator1 begin1 () {
1.571 + return find1 (0, 0, 0);
1.572 + }
1.573 + BOOST_UBLAS_INLINE
1.574 + iterator1 end1 () {
1.575 + return find1 (0, size1_, 0);
1.576 + }
1.577 +
1.578 +#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1.579 + class const_iterator2:
1.580 + public container_const_reference<triangular_matrix>,
1.581 + public random_access_iterator_base<packed_random_access_iterator_tag,
1.582 + const_iterator2, value_type> {
1.583 + public:
1.584 + typedef typename triangular_matrix::value_type value_type;
1.585 + typedef typename triangular_matrix::difference_type difference_type;
1.586 + typedef typename triangular_matrix::const_reference reference;
1.587 + typedef const typename triangular_matrix::pointer pointer;
1.588 +
1.589 + typedef const_iterator1 dual_iterator_type;
1.590 + typedef const_reverse_iterator1 dual_reverse_iterator_type;
1.591 +
1.592 + // Construction and destruction
1.593 + BOOST_UBLAS_INLINE
1.594 + const_iterator2 ():
1.595 + container_const_reference<self_type> (), it1_ (), it2_ () {}
1.596 + BOOST_UBLAS_INLINE
1.597 + const_iterator2 (const self_type &m, size_type it1, size_type it2):
1.598 + container_const_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
1.599 + BOOST_UBLAS_INLINE
1.600 + const_iterator2 (const iterator2 &it):
1.601 + container_const_reference<self_type> (it ()), it1_ (it.it1_), it2_ (it.it2_) {}
1.602 +
1.603 + // Arithmetic
1.604 + BOOST_UBLAS_INLINE
1.605 + const_iterator2 &operator ++ () {
1.606 + ++ it2_;
1.607 + return *this;
1.608 + }
1.609 + BOOST_UBLAS_INLINE
1.610 + const_iterator2 &operator -- () {
1.611 + -- it2_;
1.612 + return *this;
1.613 + }
1.614 + BOOST_UBLAS_INLINE
1.615 + const_iterator2 &operator += (difference_type n) {
1.616 + it2_ += n;
1.617 + return *this;
1.618 + }
1.619 + BOOST_UBLAS_INLINE
1.620 + const_iterator2 &operator -= (difference_type n) {
1.621 + it2_ -= n;
1.622 + return *this;
1.623 + }
1.624 + BOOST_UBLAS_INLINE
1.625 + difference_type operator - (const const_iterator2 &it) const {
1.626 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.627 + BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
1.628 + return it2_ - it.it2_;
1.629 + }
1.630 +
1.631 + // Dereference
1.632 + BOOST_UBLAS_INLINE
1.633 + const_reference operator * () const {
1.634 + return (*this) () (it1_, it2_);
1.635 + }
1.636 + BOOST_UBLAS_INLINE
1.637 + const_reference operator [] (difference_type n) const {
1.638 + return *(*this + n);
1.639 + }
1.640 +
1.641 +#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1.642 + BOOST_UBLAS_INLINE
1.643 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.644 + typename self_type::
1.645 +#endif
1.646 + const_iterator1 begin () const {
1.647 + return (*this) ().find1 (1, 0, it2_);
1.648 + }
1.649 + BOOST_UBLAS_INLINE
1.650 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.651 + typename self_type::
1.652 +#endif
1.653 + const_iterator1 end () const {
1.654 + return (*this) ().find1 (1, (*this) ().size1 (), it2_);
1.655 + }
1.656 + BOOST_UBLAS_INLINE
1.657 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.658 + typename self_type::
1.659 +#endif
1.660 + const_reverse_iterator1 rbegin () const {
1.661 + return const_reverse_iterator1 (end ());
1.662 + }
1.663 + BOOST_UBLAS_INLINE
1.664 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.665 + typename self_type::
1.666 +#endif
1.667 + const_reverse_iterator1 rend () const {
1.668 + return const_reverse_iterator1 (begin ());
1.669 + }
1.670 +#endif
1.671 +
1.672 + // Indices
1.673 + BOOST_UBLAS_INLINE
1.674 + size_type index1 () const {
1.675 + return it1_;
1.676 + }
1.677 + BOOST_UBLAS_INLINE
1.678 + size_type index2 () const {
1.679 + return it2_;
1.680 + }
1.681 +
1.682 + // Assignment
1.683 + BOOST_UBLAS_INLINE
1.684 + const_iterator2 &operator = (const const_iterator2 &it) {
1.685 + container_const_reference<self_type>::assign (&it ());
1.686 + it1_ = it.it1_;
1.687 + it2_ = it.it2_;
1.688 + return *this;
1.689 + }
1.690 +
1.691 + // Comparison
1.692 + BOOST_UBLAS_INLINE
1.693 + bool operator == (const const_iterator2 &it) const {
1.694 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.695 + BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
1.696 + return it2_ == it.it2_;
1.697 + }
1.698 + BOOST_UBLAS_INLINE
1.699 + bool operator < (const const_iterator2 &it) const {
1.700 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.701 + BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
1.702 + return it2_ < it.it2_;
1.703 + }
1.704 +
1.705 + private:
1.706 + size_type it1_;
1.707 + size_type it2_;
1.708 + };
1.709 +#endif
1.710 +
1.711 + BOOST_UBLAS_INLINE
1.712 + const_iterator2 begin2 () const {
1.713 + return find2 (0, 0, 0);
1.714 + }
1.715 + BOOST_UBLAS_INLINE
1.716 + const_iterator2 end2 () const {
1.717 + return find2 (0, 0, size2_);
1.718 + }
1.719 +
1.720 +#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1.721 + class iterator2:
1.722 + public container_reference<triangular_matrix>,
1.723 + public random_access_iterator_base<packed_random_access_iterator_tag,
1.724 + iterator2, value_type> {
1.725 + public:
1.726 + typedef typename triangular_matrix::value_type value_type;
1.727 + typedef typename triangular_matrix::difference_type difference_type;
1.728 + typedef typename triangular_matrix::reference reference;
1.729 + typedef typename triangular_matrix::pointer pointer;
1.730 +
1.731 + typedef iterator1 dual_iterator_type;
1.732 + typedef reverse_iterator1 dual_reverse_iterator_type;
1.733 +
1.734 + // Construction and destruction
1.735 + BOOST_UBLAS_INLINE
1.736 + iterator2 ():
1.737 + container_reference<self_type> (), it1_ (), it2_ () {}
1.738 + BOOST_UBLAS_INLINE
1.739 + iterator2 (self_type &m, size_type it1, size_type it2):
1.740 + container_reference<self_type> (m), it1_ (it1), it2_ (it2) {}
1.741 +
1.742 + // Arithmetic
1.743 + BOOST_UBLAS_INLINE
1.744 + iterator2 &operator ++ () {
1.745 + ++ it2_;
1.746 + return *this;
1.747 + }
1.748 + BOOST_UBLAS_INLINE
1.749 + iterator2 &operator -- () {
1.750 + -- it2_;
1.751 + return *this;
1.752 + }
1.753 + BOOST_UBLAS_INLINE
1.754 + iterator2 &operator += (difference_type n) {
1.755 + it2_ += n;
1.756 + return *this;
1.757 + }
1.758 + BOOST_UBLAS_INLINE
1.759 + iterator2 &operator -= (difference_type n) {
1.760 + it2_ -= n;
1.761 + return *this;
1.762 + }
1.763 + BOOST_UBLAS_INLINE
1.764 + difference_type operator - (const iterator2 &it) const {
1.765 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.766 + BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
1.767 + return it2_ - it.it2_;
1.768 + }
1.769 +
1.770 + // Dereference
1.771 + BOOST_UBLAS_INLINE
1.772 + reference operator * () const {
1.773 + return (*this) () (it1_, it2_);
1.774 + }
1.775 + BOOST_UBLAS_INLINE
1.776 + reference operator [] (difference_type n) const {
1.777 + return *(*this + n);
1.778 + }
1.779 +
1.780 +#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1.781 + BOOST_UBLAS_INLINE
1.782 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.783 + typename self_type::
1.784 +#endif
1.785 + iterator1 begin () const {
1.786 + return (*this) ().find1 (1, 0, it2_);
1.787 + }
1.788 + BOOST_UBLAS_INLINE
1.789 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.790 + typename self_type::
1.791 +#endif
1.792 + iterator1 end () const {
1.793 + return (*this) ().find1 (1, (*this) ().size1 (), it2_);
1.794 + }
1.795 + BOOST_UBLAS_INLINE
1.796 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.797 + typename self_type::
1.798 +#endif
1.799 + reverse_iterator1 rbegin () const {
1.800 + return reverse_iterator1 (end ());
1.801 + }
1.802 + BOOST_UBLAS_INLINE
1.803 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.804 + typename self_type::
1.805 +#endif
1.806 + reverse_iterator1 rend () const {
1.807 + return reverse_iterator1 (begin ());
1.808 + }
1.809 +#endif
1.810 +
1.811 + // Indices
1.812 + BOOST_UBLAS_INLINE
1.813 + size_type index1 () const {
1.814 + return it1_;
1.815 + }
1.816 + BOOST_UBLAS_INLINE
1.817 + size_type index2 () const {
1.818 + return it2_;
1.819 + }
1.820 +
1.821 + // Assignment
1.822 + BOOST_UBLAS_INLINE
1.823 + iterator2 &operator = (const iterator2 &it) {
1.824 + container_reference<self_type>::assign (&it ());
1.825 + it1_ = it.it1_;
1.826 + it2_ = it.it2_;
1.827 + return *this;
1.828 + }
1.829 +
1.830 + // Comparison
1.831 + BOOST_UBLAS_INLINE
1.832 + bool operator == (const iterator2 &it) const {
1.833 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.834 + BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
1.835 + return it2_ == it.it2_;
1.836 + }
1.837 + BOOST_UBLAS_INLINE
1.838 + bool operator < (const iterator2 &it) const {
1.839 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.840 + BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ());
1.841 + return it2_ < it.it2_;
1.842 + }
1.843 +
1.844 + private:
1.845 + size_type it1_;
1.846 + size_type it2_;
1.847 +
1.848 + friend class const_iterator2;
1.849 + };
1.850 +#endif
1.851 +
1.852 + BOOST_UBLAS_INLINE
1.853 + iterator2 begin2 () {
1.854 + return find2 (0, 0, 0);
1.855 + }
1.856 + BOOST_UBLAS_INLINE
1.857 + iterator2 end2 () {
1.858 + return find2 (0, 0, size2_);
1.859 + }
1.860 +
1.861 + // Reverse iterators
1.862 +
1.863 + BOOST_UBLAS_INLINE
1.864 + const_reverse_iterator1 rbegin1 () const {
1.865 + return const_reverse_iterator1 (end1 ());
1.866 + }
1.867 + BOOST_UBLAS_INLINE
1.868 + const_reverse_iterator1 rend1 () const {
1.869 + return const_reverse_iterator1 (begin1 ());
1.870 + }
1.871 +
1.872 + BOOST_UBLAS_INLINE
1.873 + reverse_iterator1 rbegin1 () {
1.874 + return reverse_iterator1 (end1 ());
1.875 + }
1.876 + BOOST_UBLAS_INLINE
1.877 + reverse_iterator1 rend1 () {
1.878 + return reverse_iterator1 (begin1 ());
1.879 + }
1.880 +
1.881 + BOOST_UBLAS_INLINE
1.882 + const_reverse_iterator2 rbegin2 () const {
1.883 + return const_reverse_iterator2 (end2 ());
1.884 + }
1.885 + BOOST_UBLAS_INLINE
1.886 + const_reverse_iterator2 rend2 () const {
1.887 + return const_reverse_iterator2 (begin2 ());
1.888 + }
1.889 +
1.890 + BOOST_UBLAS_INLINE
1.891 + reverse_iterator2 rbegin2 () {
1.892 + return reverse_iterator2 (end2 ());
1.893 + }
1.894 + BOOST_UBLAS_INLINE
1.895 + reverse_iterator2 rend2 () {
1.896 + return reverse_iterator2 (begin2 ());
1.897 + }
1.898 +
1.899 + private:
1.900 + size_type size1_;
1.901 + size_type size2_;
1.902 + array_type data_;
1.903 + static const value_type zero_;
1.904 + static const value_type one_;
1.905 + };
1.906 +
1.907 + template<class T, class TRI, class L, class A>
1.908 + const typename triangular_matrix<T, TRI, L, A>::value_type triangular_matrix<T, TRI, L, A>::zero_ = value_type/*zero*/();
1.909 + template<class T, class TRI, class L, class A>
1.910 + const typename triangular_matrix<T, TRI, L, A>::value_type triangular_matrix<T, TRI, L, A>::one_ (1);
1.911 +
1.912 +
1.913 + // Triangular matrix adaptor class
1.914 + template<class M, class TRI>
1.915 + class triangular_adaptor:
1.916 + public matrix_expression<triangular_adaptor<M, TRI> > {
1.917 +
1.918 + typedef triangular_adaptor<M, TRI> self_type;
1.919 +
1.920 + public:
1.921 +#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
1.922 + using matrix_expression<self_type>::operator ();
1.923 +#endif
1.924 + typedef const M const_matrix_type;
1.925 + typedef M matrix_type;
1.926 + typedef TRI triangular_type;
1.927 + typedef typename M::size_type size_type;
1.928 + typedef typename M::difference_type difference_type;
1.929 + typedef typename M::value_type value_type;
1.930 + typedef typename M::const_reference const_reference;
1.931 + typedef typename boost::mpl::if_<boost::is_const<M>,
1.932 + typename M::const_reference,
1.933 + typename M::reference>::type reference;
1.934 + typedef typename boost::mpl::if_<boost::is_const<M>,
1.935 + typename M::const_closure_type,
1.936 + typename M::closure_type>::type matrix_closure_type;
1.937 + typedef const self_type const_closure_type;
1.938 + typedef self_type closure_type;
1.939 + // Replaced by _temporary_traits to avoid type requirements on M
1.940 + //typedef typename M::vector_temporary_type vector_temporary_type;
1.941 + //typedef typename M::matrix_temporary_type matrix_temporary_type;
1.942 + typedef typename storage_restrict_traits<typename M::storage_category,
1.943 + packed_proxy_tag>::storage_category storage_category;
1.944 + typedef typename M::orientation_category orientation_category;
1.945 +
1.946 + // Construction and destruction
1.947 + BOOST_UBLAS_INLINE
1.948 + triangular_adaptor (matrix_type &data):
1.949 + matrix_expression<self_type> (),
1.950 + data_ (data) {}
1.951 + BOOST_UBLAS_INLINE
1.952 + triangular_adaptor (const triangular_adaptor &m):
1.953 + matrix_expression<self_type> (),
1.954 + data_ (m.data_) {}
1.955 +
1.956 + // Accessors
1.957 + BOOST_UBLAS_INLINE
1.958 + size_type size1 () const {
1.959 + return data_.size1 ();
1.960 + }
1.961 + BOOST_UBLAS_INLINE
1.962 + size_type size2 () const {
1.963 + return data_.size2 ();
1.964 + }
1.965 +
1.966 + // Storage accessors
1.967 + BOOST_UBLAS_INLINE
1.968 + const matrix_closure_type &data () const {
1.969 + return data_;
1.970 + }
1.971 + BOOST_UBLAS_INLINE
1.972 + matrix_closure_type &data () {
1.973 + return data_;
1.974 + }
1.975 +
1.976 + // Element access
1.977 +#ifndef BOOST_UBLAS_PROXY_CONST_MEMBER
1.978 + BOOST_UBLAS_INLINE
1.979 + const_reference operator () (size_type i, size_type j) const {
1.980 + BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
1.981 + BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
1.982 + if (triangular_type::other (i, j))
1.983 + return data () (i, j);
1.984 + else if (triangular_type::one (i, j))
1.985 + return one_;
1.986 + else
1.987 + return zero_;
1.988 + }
1.989 + BOOST_UBLAS_INLINE
1.990 + reference operator () (size_type i, size_type j) {
1.991 + BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
1.992 + BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
1.993 + if (triangular_type::other (i, j))
1.994 + return data () (i, j);
1.995 + else if (triangular_type::one (i, j)) {
1.996 + bad_index ().raise ();
1.997 + // arbitary return value
1.998 + return const_cast<reference>(one_);
1.999 + } else {
1.1000 + bad_index ().raise ();
1.1001 + // arbitary return value
1.1002 + return const_cast<reference>(zero_);
1.1003 + }
1.1004 + }
1.1005 +#else
1.1006 + BOOST_UBLAS_INLINE
1.1007 + reference operator () (size_type i, size_type j) const {
1.1008 + BOOST_UBLAS_CHECK (i < size1 (), bad_index ());
1.1009 + BOOST_UBLAS_CHECK (j < size2 (), bad_index ());
1.1010 + if (triangular_type::other (i, j))
1.1011 + return data () (i, j);
1.1012 + else if (triangular_type::one (i, j)) {
1.1013 + bad_index ().raise ();
1.1014 + // arbitary return value
1.1015 + return const_cast<reference>(one_);
1.1016 + } else {
1.1017 + bad_index ().raise ();
1.1018 + // arbitary return value
1.1019 + return const_cast<reference>(zero_);
1.1020 + }
1.1021 + }
1.1022 +#endif
1.1023 +
1.1024 + // Assignment
1.1025 + BOOST_UBLAS_INLINE
1.1026 + triangular_adaptor &operator = (const triangular_adaptor &m) {
1.1027 + matrix_assign<scalar_assign> (*this, m);
1.1028 + return *this;
1.1029 + }
1.1030 + BOOST_UBLAS_INLINE
1.1031 + triangular_adaptor &assign_temporary (triangular_adaptor &m) {
1.1032 + *this = m;
1.1033 + return *this;
1.1034 + }
1.1035 + template<class AE>
1.1036 + BOOST_UBLAS_INLINE
1.1037 + triangular_adaptor &operator = (const matrix_expression<AE> &ae) {
1.1038 + matrix_assign<scalar_assign> (*this, matrix<value_type> (ae));
1.1039 + return *this;
1.1040 + }
1.1041 + template<class AE>
1.1042 + BOOST_UBLAS_INLINE
1.1043 + triangular_adaptor &assign (const matrix_expression<AE> &ae) {
1.1044 + matrix_assign<scalar_assign> (*this, ae);
1.1045 + return *this;
1.1046 + }
1.1047 + template<class AE>
1.1048 + BOOST_UBLAS_INLINE
1.1049 + triangular_adaptor& operator += (const matrix_expression<AE> &ae) {
1.1050 + matrix_assign<scalar_assign> (*this, matrix<value_type> (*this + ae));
1.1051 + return *this;
1.1052 + }
1.1053 + template<class AE>
1.1054 + BOOST_UBLAS_INLINE
1.1055 + triangular_adaptor &plus_assign (const matrix_expression<AE> &ae) {
1.1056 + matrix_assign<scalar_plus_assign> (*this, ae);
1.1057 + return *this;
1.1058 + }
1.1059 + template<class AE>
1.1060 + BOOST_UBLAS_INLINE
1.1061 + triangular_adaptor& operator -= (const matrix_expression<AE> &ae) {
1.1062 + matrix_assign<scalar_assign> (*this, matrix<value_type> (*this - ae));
1.1063 + return *this;
1.1064 + }
1.1065 + template<class AE>
1.1066 + BOOST_UBLAS_INLINE
1.1067 + triangular_adaptor &minus_assign (const matrix_expression<AE> &ae) {
1.1068 + matrix_assign<scalar_minus_assign> (*this, ae);
1.1069 + return *this;
1.1070 + }
1.1071 + template<class AT>
1.1072 + BOOST_UBLAS_INLINE
1.1073 + triangular_adaptor& operator *= (const AT &at) {
1.1074 + matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
1.1075 + return *this;
1.1076 + }
1.1077 + template<class AT>
1.1078 + BOOST_UBLAS_INLINE
1.1079 + triangular_adaptor& operator /= (const AT &at) {
1.1080 + matrix_assign_scalar<scalar_divides_assign> (*this, at);
1.1081 + return *this;
1.1082 + }
1.1083 +
1.1084 + // Closure comparison
1.1085 + BOOST_UBLAS_INLINE
1.1086 + bool same_closure (const triangular_adaptor &ta) const {
1.1087 + return (*this).data ().same_closure (ta.data ());
1.1088 + }
1.1089 +
1.1090 + // Swapping
1.1091 + BOOST_UBLAS_INLINE
1.1092 + void swap (triangular_adaptor &m) {
1.1093 + if (this != &m)
1.1094 + matrix_swap<scalar_swap> (*this, m);
1.1095 + }
1.1096 + BOOST_UBLAS_INLINE
1.1097 + friend void swap (triangular_adaptor &m1, triangular_adaptor &m2) {
1.1098 + m1.swap (m2);
1.1099 + }
1.1100 +
1.1101 + // Iterator types
1.1102 + private:
1.1103 + typedef typename M::const_iterator1 const_subiterator1_type;
1.1104 + typedef typename boost::mpl::if_<boost::is_const<M>,
1.1105 + typename M::const_iterator1,
1.1106 + typename M::iterator1>::type subiterator1_type;
1.1107 + typedef typename M::const_iterator2 const_subiterator2_type;
1.1108 + typedef typename boost::mpl::if_<boost::is_const<M>,
1.1109 + typename M::const_iterator2,
1.1110 + typename M::iterator2>::type subiterator2_type;
1.1111 +
1.1112 + public:
1.1113 +#ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR
1.1114 + typedef indexed_iterator1<self_type, packed_random_access_iterator_tag> iterator1;
1.1115 + typedef indexed_iterator2<self_type, packed_random_access_iterator_tag> iterator2;
1.1116 + typedef indexed_const_iterator1<self_type, packed_random_access_iterator_tag> const_iterator1;
1.1117 + typedef indexed_const_iterator2<self_type, packed_random_access_iterator_tag> const_iterator2;
1.1118 +#else
1.1119 + class const_iterator1;
1.1120 + class iterator1;
1.1121 + class const_iterator2;
1.1122 + class iterator2;
1.1123 +#endif
1.1124 + typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1;
1.1125 + typedef reverse_iterator_base1<iterator1> reverse_iterator1;
1.1126 + typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2;
1.1127 + typedef reverse_iterator_base2<iterator2> reverse_iterator2;
1.1128 +
1.1129 + // Element lookup
1.1130 + BOOST_UBLAS_INLINE
1.1131 + const_iterator1 find1 (int rank, size_type i, size_type j) const {
1.1132 + if (rank == 1)
1.1133 + i = triangular_type::restrict1 (i, j);
1.1134 + return const_iterator1 (*this, data ().find1 (rank, i, j));
1.1135 + }
1.1136 + BOOST_UBLAS_INLINE
1.1137 + iterator1 find1 (int rank, size_type i, size_type j) {
1.1138 + if (rank == 1)
1.1139 + i = triangular_type::mutable_restrict1 (i, j);
1.1140 + return iterator1 (*this, data ().find1 (rank, i, j));
1.1141 + }
1.1142 + BOOST_UBLAS_INLINE
1.1143 + const_iterator2 find2 (int rank, size_type i, size_type j) const {
1.1144 + if (rank == 1)
1.1145 + j = triangular_type::restrict2 (i, j);
1.1146 + return const_iterator2 (*this, data ().find2 (rank, i, j));
1.1147 + }
1.1148 + BOOST_UBLAS_INLINE
1.1149 + iterator2 find2 (int rank, size_type i, size_type j) {
1.1150 + if (rank == 1)
1.1151 + j = triangular_type::mutable_restrict2 (i, j);
1.1152 + return iterator2 (*this, data ().find2 (rank, i, j));
1.1153 + }
1.1154 +
1.1155 + // Iterators simply are indices.
1.1156 +
1.1157 +#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1.1158 + class const_iterator1:
1.1159 + public container_const_reference<triangular_adaptor>,
1.1160 + public random_access_iterator_base<typename iterator_restrict_traits<
1.1161 + typename const_subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,
1.1162 + const_iterator1, value_type> {
1.1163 + public:
1.1164 + typedef typename const_subiterator1_type::value_type value_type;
1.1165 + typedef typename const_subiterator1_type::difference_type difference_type;
1.1166 + typedef typename const_subiterator1_type::reference reference;
1.1167 + typedef typename const_subiterator1_type::pointer pointer;
1.1168 +
1.1169 + typedef const_iterator2 dual_iterator_type;
1.1170 + typedef const_reverse_iterator2 dual_reverse_iterator_type;
1.1171 +
1.1172 + // Construction and destruction
1.1173 + BOOST_UBLAS_INLINE
1.1174 + const_iterator1 ():
1.1175 + container_const_reference<self_type> (), it1_ () {}
1.1176 + BOOST_UBLAS_INLINE
1.1177 + const_iterator1 (const self_type &m, const const_subiterator1_type &it1):
1.1178 + container_const_reference<self_type> (m), it1_ (it1) {}
1.1179 + BOOST_UBLAS_INLINE
1.1180 + const_iterator1 (const iterator1 &it):
1.1181 + container_const_reference<self_type> (it ()), it1_ (it.it1_) {}
1.1182 +
1.1183 + // Arithmetic
1.1184 + BOOST_UBLAS_INLINE
1.1185 + const_iterator1 &operator ++ () {
1.1186 + ++ it1_;
1.1187 + return *this;
1.1188 + }
1.1189 + BOOST_UBLAS_INLINE
1.1190 + const_iterator1 &operator -- () {
1.1191 + -- it1_;
1.1192 + return *this;
1.1193 + }
1.1194 + BOOST_UBLAS_INLINE
1.1195 + const_iterator1 &operator += (difference_type n) {
1.1196 + it1_ += n;
1.1197 + return *this;
1.1198 + }
1.1199 + BOOST_UBLAS_INLINE
1.1200 + const_iterator1 &operator -= (difference_type n) {
1.1201 + it1_ -= n;
1.1202 + return *this;
1.1203 + }
1.1204 + BOOST_UBLAS_INLINE
1.1205 + difference_type operator - (const const_iterator1 &it) const {
1.1206 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.1207 + return it1_ - it.it1_;
1.1208 + }
1.1209 +
1.1210 + // Dereference
1.1211 + BOOST_UBLAS_INLINE
1.1212 + const_reference operator * () const {
1.1213 + size_type i = index1 ();
1.1214 + size_type j = index2 ();
1.1215 + BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());
1.1216 + BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());
1.1217 + if (triangular_type::other (i, j))
1.1218 + return *it1_;
1.1219 + else
1.1220 + return (*this) () (i, j);
1.1221 + }
1.1222 + BOOST_UBLAS_INLINE
1.1223 + const_reference operator [] (difference_type n) const {
1.1224 + return *(*this + n);
1.1225 + }
1.1226 +
1.1227 +#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1.1228 + BOOST_UBLAS_INLINE
1.1229 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.1230 + typename self_type::
1.1231 +#endif
1.1232 + const_iterator2 begin () const {
1.1233 + return (*this) ().find2 (1, index1 (), 0);
1.1234 + }
1.1235 + BOOST_UBLAS_INLINE
1.1236 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.1237 + typename self_type::
1.1238 +#endif
1.1239 + const_iterator2 end () const {
1.1240 + return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
1.1241 + }
1.1242 + BOOST_UBLAS_INLINE
1.1243 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.1244 + typename self_type::
1.1245 +#endif
1.1246 + const_reverse_iterator2 rbegin () const {
1.1247 + return const_reverse_iterator2 (end ());
1.1248 + }
1.1249 + BOOST_UBLAS_INLINE
1.1250 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.1251 + typename self_type::
1.1252 +#endif
1.1253 + const_reverse_iterator2 rend () const {
1.1254 + return const_reverse_iterator2 (begin ());
1.1255 + }
1.1256 +#endif
1.1257 +
1.1258 + // Indices
1.1259 + BOOST_UBLAS_INLINE
1.1260 + size_type index1 () const {
1.1261 + return it1_.index1 ();
1.1262 + }
1.1263 + BOOST_UBLAS_INLINE
1.1264 + size_type index2 () const {
1.1265 + return it1_.index2 ();
1.1266 + }
1.1267 +
1.1268 + // Assignment
1.1269 + BOOST_UBLAS_INLINE
1.1270 + const_iterator1 &operator = (const const_iterator1 &it) {
1.1271 + container_const_reference<self_type>::assign (&it ());
1.1272 + it1_ = it.it1_;
1.1273 + return *this;
1.1274 + }
1.1275 +
1.1276 + // Comparison
1.1277 + BOOST_UBLAS_INLINE
1.1278 + bool operator == (const const_iterator1 &it) const {
1.1279 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.1280 + return it1_ == it.it1_;
1.1281 + }
1.1282 + BOOST_UBLAS_INLINE
1.1283 + bool operator < (const const_iterator1 &it) const {
1.1284 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.1285 + return it1_ < it.it1_;
1.1286 + }
1.1287 +
1.1288 + private:
1.1289 + const_subiterator1_type it1_;
1.1290 + };
1.1291 +#endif
1.1292 +
1.1293 + BOOST_UBLAS_INLINE
1.1294 + const_iterator1 begin1 () const {
1.1295 + return find1 (0, 0, 0);
1.1296 + }
1.1297 + BOOST_UBLAS_INLINE
1.1298 + const_iterator1 end1 () const {
1.1299 + return find1 (0, size1 (), 0);
1.1300 + }
1.1301 +
1.1302 +#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1.1303 + class iterator1:
1.1304 + public container_reference<triangular_adaptor>,
1.1305 + public random_access_iterator_base<typename iterator_restrict_traits<
1.1306 + typename subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,
1.1307 + iterator1, value_type> {
1.1308 + public:
1.1309 + typedef typename subiterator1_type::value_type value_type;
1.1310 + typedef typename subiterator1_type::difference_type difference_type;
1.1311 + typedef typename subiterator1_type::reference reference;
1.1312 + typedef typename subiterator1_type::pointer pointer;
1.1313 +
1.1314 + typedef iterator2 dual_iterator_type;
1.1315 + typedef reverse_iterator2 dual_reverse_iterator_type;
1.1316 +
1.1317 + // Construction and destruction
1.1318 + BOOST_UBLAS_INLINE
1.1319 + iterator1 ():
1.1320 + container_reference<self_type> (), it1_ () {}
1.1321 + BOOST_UBLAS_INLINE
1.1322 + iterator1 (self_type &m, const subiterator1_type &it1):
1.1323 + container_reference<self_type> (m), it1_ (it1) {}
1.1324 +
1.1325 + // Arithmetic
1.1326 + BOOST_UBLAS_INLINE
1.1327 + iterator1 &operator ++ () {
1.1328 + ++ it1_;
1.1329 + return *this;
1.1330 + }
1.1331 + BOOST_UBLAS_INLINE
1.1332 + iterator1 &operator -- () {
1.1333 + -- it1_;
1.1334 + return *this;
1.1335 + }
1.1336 + BOOST_UBLAS_INLINE
1.1337 + iterator1 &operator += (difference_type n) {
1.1338 + it1_ += n;
1.1339 + return *this;
1.1340 + }
1.1341 + BOOST_UBLAS_INLINE
1.1342 + iterator1 &operator -= (difference_type n) {
1.1343 + it1_ -= n;
1.1344 + return *this;
1.1345 + }
1.1346 + BOOST_UBLAS_INLINE
1.1347 + difference_type operator - (const iterator1 &it) const {
1.1348 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.1349 + return it1_ - it.it1_;
1.1350 + }
1.1351 +
1.1352 + // Dereference
1.1353 + BOOST_UBLAS_INLINE
1.1354 + reference operator * () const {
1.1355 + size_type i = index1 ();
1.1356 + size_type j = index2 ();
1.1357 + BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());
1.1358 + BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());
1.1359 + if (triangular_type::other (i, j))
1.1360 + return *it1_;
1.1361 + else
1.1362 + return (*this) () (i, j);
1.1363 + }
1.1364 + BOOST_UBLAS_INLINE
1.1365 + reference operator [] (difference_type n) const {
1.1366 + return *(*this + n);
1.1367 + }
1.1368 +
1.1369 +#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1.1370 + BOOST_UBLAS_INLINE
1.1371 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.1372 + typename self_type::
1.1373 +#endif
1.1374 + iterator2 begin () const {
1.1375 + return (*this) ().find2 (1, index1 (), 0);
1.1376 + }
1.1377 + BOOST_UBLAS_INLINE
1.1378 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.1379 + typename self_type::
1.1380 +#endif
1.1381 + iterator2 end () const {
1.1382 + return (*this) ().find2 (1, index1 (), (*this) ().size2 ());
1.1383 + }
1.1384 + BOOST_UBLAS_INLINE
1.1385 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.1386 + typename self_type::
1.1387 +#endif
1.1388 + reverse_iterator2 rbegin () const {
1.1389 + return reverse_iterator2 (end ());
1.1390 + }
1.1391 + BOOST_UBLAS_INLINE
1.1392 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.1393 + typename self_type::
1.1394 +#endif
1.1395 + reverse_iterator2 rend () const {
1.1396 + return reverse_iterator2 (begin ());
1.1397 + }
1.1398 +#endif
1.1399 +
1.1400 + // Indices
1.1401 + BOOST_UBLAS_INLINE
1.1402 + size_type index1 () const {
1.1403 + return it1_.index1 ();
1.1404 + }
1.1405 + BOOST_UBLAS_INLINE
1.1406 + size_type index2 () const {
1.1407 + return it1_.index2 ();
1.1408 + }
1.1409 +
1.1410 + // Assignment
1.1411 + BOOST_UBLAS_INLINE
1.1412 + iterator1 &operator = (const iterator1 &it) {
1.1413 + container_reference<self_type>::assign (&it ());
1.1414 + it1_ = it.it1_;
1.1415 + return *this;
1.1416 + }
1.1417 +
1.1418 + // Comparison
1.1419 + BOOST_UBLAS_INLINE
1.1420 + bool operator == (const iterator1 &it) const {
1.1421 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.1422 + return it1_ == it.it1_;
1.1423 + }
1.1424 + BOOST_UBLAS_INLINE
1.1425 + bool operator < (const iterator1 &it) const {
1.1426 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.1427 + return it1_ < it.it1_;
1.1428 + }
1.1429 +
1.1430 + private:
1.1431 + subiterator1_type it1_;
1.1432 +
1.1433 + friend class const_iterator1;
1.1434 + };
1.1435 +#endif
1.1436 +
1.1437 + BOOST_UBLAS_INLINE
1.1438 + iterator1 begin1 () {
1.1439 + return find1 (0, 0, 0);
1.1440 + }
1.1441 + BOOST_UBLAS_INLINE
1.1442 + iterator1 end1 () {
1.1443 + return find1 (0, size1 (), 0);
1.1444 + }
1.1445 +
1.1446 +#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1.1447 + class const_iterator2:
1.1448 + public container_const_reference<triangular_adaptor>,
1.1449 + public random_access_iterator_base<typename iterator_restrict_traits<
1.1450 + typename const_subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,
1.1451 + const_iterator2, value_type> {
1.1452 + public:
1.1453 + typedef typename const_subiterator2_type::value_type value_type;
1.1454 + typedef typename const_subiterator2_type::difference_type difference_type;
1.1455 + typedef typename const_subiterator2_type::reference reference;
1.1456 + typedef typename const_subiterator2_type::pointer pointer;
1.1457 +
1.1458 + typedef const_iterator1 dual_iterator_type;
1.1459 + typedef const_reverse_iterator1 dual_reverse_iterator_type;
1.1460 +
1.1461 + // Construction and destruction
1.1462 + BOOST_UBLAS_INLINE
1.1463 + const_iterator2 ():
1.1464 + container_const_reference<self_type> (), it2_ () {}
1.1465 + BOOST_UBLAS_INLINE
1.1466 + const_iterator2 (const self_type &m, const const_subiterator2_type &it2):
1.1467 + container_const_reference<self_type> (m), it2_ (it2) {}
1.1468 + BOOST_UBLAS_INLINE
1.1469 + const_iterator2 (const iterator2 &it):
1.1470 + container_const_reference<self_type> (it ()), it2_ (it.it2_) {}
1.1471 +
1.1472 + // Arithmetic
1.1473 + BOOST_UBLAS_INLINE
1.1474 + const_iterator2 &operator ++ () {
1.1475 + ++ it2_;
1.1476 + return *this;
1.1477 + }
1.1478 + BOOST_UBLAS_INLINE
1.1479 + const_iterator2 &operator -- () {
1.1480 + -- it2_;
1.1481 + return *this;
1.1482 + }
1.1483 + BOOST_UBLAS_INLINE
1.1484 + const_iterator2 &operator += (difference_type n) {
1.1485 + it2_ += n;
1.1486 + return *this;
1.1487 + }
1.1488 + BOOST_UBLAS_INLINE
1.1489 + const_iterator2 &operator -= (difference_type n) {
1.1490 + it2_ -= n;
1.1491 + return *this;
1.1492 + }
1.1493 + BOOST_UBLAS_INLINE
1.1494 + difference_type operator - (const const_iterator2 &it) const {
1.1495 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.1496 + return it2_ - it.it2_;
1.1497 + }
1.1498 +
1.1499 + // Dereference
1.1500 + BOOST_UBLAS_INLINE
1.1501 + const_reference operator * () const {
1.1502 + size_type i = index1 ();
1.1503 + size_type j = index2 ();
1.1504 + BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());
1.1505 + BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());
1.1506 + if (triangular_type::other (i, j))
1.1507 + return *it2_;
1.1508 + else
1.1509 + return (*this) () (i, j);
1.1510 + }
1.1511 + BOOST_UBLAS_INLINE
1.1512 + const_reference operator [] (difference_type n) const {
1.1513 + return *(*this + n);
1.1514 + }
1.1515 +
1.1516 +#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1.1517 + BOOST_UBLAS_INLINE
1.1518 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.1519 + typename self_type::
1.1520 +#endif
1.1521 + const_iterator1 begin () const {
1.1522 + return (*this) ().find1 (1, 0, index2 ());
1.1523 + }
1.1524 + BOOST_UBLAS_INLINE
1.1525 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.1526 + typename self_type::
1.1527 +#endif
1.1528 + const_iterator1 end () const {
1.1529 + return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
1.1530 + }
1.1531 + BOOST_UBLAS_INLINE
1.1532 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.1533 + typename self_type::
1.1534 +#endif
1.1535 + const_reverse_iterator1 rbegin () const {
1.1536 + return const_reverse_iterator1 (end ());
1.1537 + }
1.1538 + BOOST_UBLAS_INLINE
1.1539 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.1540 + typename self_type::
1.1541 +#endif
1.1542 + const_reverse_iterator1 rend () const {
1.1543 + return const_reverse_iterator1 (begin ());
1.1544 + }
1.1545 +#endif
1.1546 +
1.1547 + // Indices
1.1548 + BOOST_UBLAS_INLINE
1.1549 + size_type index1 () const {
1.1550 + return it2_.index1 ();
1.1551 + }
1.1552 + BOOST_UBLAS_INLINE
1.1553 + size_type index2 () const {
1.1554 + return it2_.index2 ();
1.1555 + }
1.1556 +
1.1557 + // Assignment
1.1558 + BOOST_UBLAS_INLINE
1.1559 + const_iterator2 &operator = (const const_iterator2 &it) {
1.1560 + container_const_reference<self_type>::assign (&it ());
1.1561 + it2_ = it.it2_;
1.1562 + return *this;
1.1563 + }
1.1564 +
1.1565 + // Comparison
1.1566 + BOOST_UBLAS_INLINE
1.1567 + bool operator == (const const_iterator2 &it) const {
1.1568 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.1569 + return it2_ == it.it2_;
1.1570 + }
1.1571 + BOOST_UBLAS_INLINE
1.1572 + bool operator < (const const_iterator2 &it) const {
1.1573 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.1574 + return it2_ < it.it2_;
1.1575 + }
1.1576 +
1.1577 + private:
1.1578 + const_subiterator2_type it2_;
1.1579 + };
1.1580 +#endif
1.1581 +
1.1582 + BOOST_UBLAS_INLINE
1.1583 + const_iterator2 begin2 () const {
1.1584 + return find2 (0, 0, 0);
1.1585 + }
1.1586 + BOOST_UBLAS_INLINE
1.1587 + const_iterator2 end2 () const {
1.1588 + return find2 (0, 0, size2 ());
1.1589 + }
1.1590 +
1.1591 +#ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR
1.1592 + class iterator2:
1.1593 + public container_reference<triangular_adaptor>,
1.1594 + public random_access_iterator_base<typename iterator_restrict_traits<
1.1595 + typename subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category,
1.1596 + iterator2, value_type> {
1.1597 + public:
1.1598 + typedef typename subiterator2_type::value_type value_type;
1.1599 + typedef typename subiterator2_type::difference_type difference_type;
1.1600 + typedef typename subiterator2_type::reference reference;
1.1601 + typedef typename subiterator2_type::pointer pointer;
1.1602 +
1.1603 + typedef iterator1 dual_iterator_type;
1.1604 + typedef reverse_iterator1 dual_reverse_iterator_type;
1.1605 +
1.1606 + // Construction and destruction
1.1607 + BOOST_UBLAS_INLINE
1.1608 + iterator2 ():
1.1609 + container_reference<self_type> (), it2_ () {}
1.1610 + BOOST_UBLAS_INLINE
1.1611 + iterator2 (self_type &m, const subiterator2_type &it2):
1.1612 + container_reference<self_type> (m), it2_ (it2) {}
1.1613 +
1.1614 + // Arithmetic
1.1615 + BOOST_UBLAS_INLINE
1.1616 + iterator2 &operator ++ () {
1.1617 + ++ it2_;
1.1618 + return *this;
1.1619 + }
1.1620 + BOOST_UBLAS_INLINE
1.1621 + iterator2 &operator -- () {
1.1622 + -- it2_;
1.1623 + return *this;
1.1624 + }
1.1625 + BOOST_UBLAS_INLINE
1.1626 + iterator2 &operator += (difference_type n) {
1.1627 + it2_ += n;
1.1628 + return *this;
1.1629 + }
1.1630 + BOOST_UBLAS_INLINE
1.1631 + iterator2 &operator -= (difference_type n) {
1.1632 + it2_ -= n;
1.1633 + return *this;
1.1634 + }
1.1635 + BOOST_UBLAS_INLINE
1.1636 + difference_type operator - (const iterator2 &it) const {
1.1637 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.1638 + return it2_ - it.it2_;
1.1639 + }
1.1640 +
1.1641 + // Dereference
1.1642 + BOOST_UBLAS_INLINE
1.1643 + reference operator * () const {
1.1644 + size_type i = index1 ();
1.1645 + size_type j = index2 ();
1.1646 + BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ());
1.1647 + BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ());
1.1648 + if (triangular_type::other (i, j))
1.1649 + return *it2_;
1.1650 + else
1.1651 + return (*this) () (i, j);
1.1652 + }
1.1653 + BOOST_UBLAS_INLINE
1.1654 + reference operator [] (difference_type n) const {
1.1655 + return *(*this + n);
1.1656 + }
1.1657 +
1.1658 +#ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1.1659 + BOOST_UBLAS_INLINE
1.1660 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.1661 + typename self_type::
1.1662 +#endif
1.1663 + iterator1 begin () const {
1.1664 + return (*this) ().find1 (1, 0, index2 ());
1.1665 + }
1.1666 + BOOST_UBLAS_INLINE
1.1667 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.1668 + typename self_type::
1.1669 +#endif
1.1670 + iterator1 end () const {
1.1671 + return (*this) ().find1 (1, (*this) ().size1 (), index2 ());
1.1672 + }
1.1673 + BOOST_UBLAS_INLINE
1.1674 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.1675 + typename self_type::
1.1676 +#endif
1.1677 + reverse_iterator1 rbegin () const {
1.1678 + return reverse_iterator1 (end ());
1.1679 + }
1.1680 + BOOST_UBLAS_INLINE
1.1681 +#ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1.1682 + typename self_type::
1.1683 +#endif
1.1684 + reverse_iterator1 rend () const {
1.1685 + return reverse_iterator1 (begin ());
1.1686 + }
1.1687 +#endif
1.1688 +
1.1689 + // Indices
1.1690 + BOOST_UBLAS_INLINE
1.1691 + size_type index1 () const {
1.1692 + return it2_.index1 ();
1.1693 + }
1.1694 + BOOST_UBLAS_INLINE
1.1695 + size_type index2 () const {
1.1696 + return it2_.index2 ();
1.1697 + }
1.1698 +
1.1699 + // Assignment
1.1700 + BOOST_UBLAS_INLINE
1.1701 + iterator2 &operator = (const iterator2 &it) {
1.1702 + container_reference<self_type>::assign (&it ());
1.1703 + it2_ = it.it2_;
1.1704 + return *this;
1.1705 + }
1.1706 +
1.1707 + // Comparison
1.1708 + BOOST_UBLAS_INLINE
1.1709 + bool operator == (const iterator2 &it) const {
1.1710 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.1711 + return it2_ == it.it2_;
1.1712 + }
1.1713 + BOOST_UBLAS_INLINE
1.1714 + bool operator < (const iterator2 &it) const {
1.1715 + BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1.1716 + return it2_ < it.it2_;
1.1717 + }
1.1718 +
1.1719 + private:
1.1720 + subiterator2_type it2_;
1.1721 +
1.1722 + friend class const_iterator2;
1.1723 + };
1.1724 +#endif
1.1725 +
1.1726 + BOOST_UBLAS_INLINE
1.1727 + iterator2 begin2 () {
1.1728 + return find2 (0, 0, 0);
1.1729 + }
1.1730 + BOOST_UBLAS_INLINE
1.1731 + iterator2 end2 () {
1.1732 + return find2 (0, 0, size2 ());
1.1733 + }
1.1734 +
1.1735 + // Reverse iterators
1.1736 +
1.1737 + BOOST_UBLAS_INLINE
1.1738 + const_reverse_iterator1 rbegin1 () const {
1.1739 + return const_reverse_iterator1 (end1 ());
1.1740 + }
1.1741 + BOOST_UBLAS_INLINE
1.1742 + const_reverse_iterator1 rend1 () const {
1.1743 + return const_reverse_iterator1 (begin1 ());
1.1744 + }
1.1745 +
1.1746 + BOOST_UBLAS_INLINE
1.1747 + reverse_iterator1 rbegin1 () {
1.1748 + return reverse_iterator1 (end1 ());
1.1749 + }
1.1750 + BOOST_UBLAS_INLINE
1.1751 + reverse_iterator1 rend1 () {
1.1752 + return reverse_iterator1 (begin1 ());
1.1753 + }
1.1754 +
1.1755 + BOOST_UBLAS_INLINE
1.1756 + const_reverse_iterator2 rbegin2 () const {
1.1757 + return const_reverse_iterator2 (end2 ());
1.1758 + }
1.1759 + BOOST_UBLAS_INLINE
1.1760 + const_reverse_iterator2 rend2 () const {
1.1761 + return const_reverse_iterator2 (begin2 ());
1.1762 + }
1.1763 +
1.1764 + BOOST_UBLAS_INLINE
1.1765 + reverse_iterator2 rbegin2 () {
1.1766 + return reverse_iterator2 (end2 ());
1.1767 + }
1.1768 + BOOST_UBLAS_INLINE
1.1769 + reverse_iterator2 rend2 () {
1.1770 + return reverse_iterator2 (begin2 ());
1.1771 + }
1.1772 +
1.1773 + private:
1.1774 + matrix_closure_type data_;
1.1775 + static const value_type zero_;
1.1776 + static const value_type one_;
1.1777 + };
1.1778 +
1.1779 + template<class M, class TRI>
1.1780 + const typename triangular_adaptor<M, TRI>::value_type triangular_adaptor<M, TRI>::zero_ = value_type/*zero*/();
1.1781 + template<class M, class TRI>
1.1782 + const typename triangular_adaptor<M, TRI>::value_type triangular_adaptor<M, TRI>::one_ (1);
1.1783 +
1.1784 + template <class M, class TRI>
1.1785 + struct vector_temporary_traits< triangular_adaptor<M, TRI> >
1.1786 + : vector_temporary_traits< typename boost::remove_const<M>::type > {} ;
1.1787 + template <class M, class TRI>
1.1788 + struct vector_temporary_traits< const triangular_adaptor<M, TRI> >
1.1789 + : vector_temporary_traits< typename boost::remove_const<M>::type > {} ;
1.1790 +
1.1791 + template <class M, class TRI>
1.1792 + struct matrix_temporary_traits< triangular_adaptor<M, TRI> >
1.1793 + : matrix_temporary_traits< typename boost::remove_const<M>::type > {};
1.1794 + template <class M, class TRI>
1.1795 + struct matrix_temporary_traits< const triangular_adaptor<M, TRI> >
1.1796 + : matrix_temporary_traits< typename boost::remove_const<M>::type > {};
1.1797 +
1.1798 +
1.1799 + template<class E1, class E2>
1.1800 + struct matrix_vector_solve_traits {
1.1801 + typedef typename promote_traits<typename E1::value_type, typename E2::value_type>::promote_type promote_type;
1.1802 + typedef vector<promote_type> result_type;
1.1803 + };
1.1804 +
1.1805 + // Operations:
1.1806 + // n * (n - 1) / 2 + n = n * (n + 1) / 2 multiplications,
1.1807 + // n * (n - 1) / 2 additions
1.1808 +
1.1809 + // Dense (proxy) case
1.1810 + template<class E1, class E2>
1.1811 + BOOST_UBLAS_INLINE
1.1812 + void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
1.1813 + lower_tag, column_major_tag, dense_proxy_tag) {
1.1814 + typedef typename E2::size_type size_type;
1.1815 + typedef typename E2::difference_type difference_type;
1.1816 + typedef typename E2::value_type value_type;
1.1817 +
1.1818 + BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
1.1819 + BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
1.1820 + size_type size = e2 ().size ();
1.1821 + for (size_type n = 0; n < size; ++ n) {
1.1822 +#ifndef BOOST_UBLAS_SINGULAR_CHECK
1.1823 + BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
1.1824 +#else
1.1825 + if (e1 () (n, n) == value_type/*zero*/())
1.1826 + singular ().raise ();
1.1827 +#endif
1.1828 + value_type t = e2 () (n) /= e1 () (n, n);
1.1829 + if (t != value_type/*zero*/()) {
1.1830 + for (size_type m = n + 1; m < size; ++ m)
1.1831 + e2 () (m) -= e1 () (m, n) * t;
1.1832 + }
1.1833 + }
1.1834 + }
1.1835 + // Packed (proxy) case
1.1836 + template<class E1, class E2>
1.1837 + BOOST_UBLAS_INLINE
1.1838 + void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
1.1839 + lower_tag, column_major_tag, packed_proxy_tag) {
1.1840 + typedef typename E2::size_type size_type;
1.1841 + typedef typename E2::difference_type difference_type;
1.1842 + typedef typename E2::value_type value_type;
1.1843 +
1.1844 + BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
1.1845 + BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
1.1846 + size_type size = e2 ().size ();
1.1847 + for (size_type n = 0; n < size; ++ n) {
1.1848 +#ifndef BOOST_UBLAS_SINGULAR_CHECK
1.1849 + BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
1.1850 +#else
1.1851 + if (e1 () (n, n) == value_type/*zero*/())
1.1852 + singular ().raise ();
1.1853 +#endif
1.1854 + value_type t = e2 () (n) /= e1 () (n, n);
1.1855 + if (t != value_type/*zero*/()) {
1.1856 + typename E1::const_iterator1 it1e1 (e1 ().find1 (1, n + 1, n));
1.1857 + typename E1::const_iterator1 it1e1_end (e1 ().find1 (1, e1 ().size1 (), n));
1.1858 + difference_type m (it1e1_end - it1e1);
1.1859 + while (-- m >= 0)
1.1860 + e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1;
1.1861 + }
1.1862 + }
1.1863 + }
1.1864 + // Sparse (proxy) case
1.1865 + template<class E1, class E2>
1.1866 + BOOST_UBLAS_INLINE
1.1867 + void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
1.1868 + lower_tag, column_major_tag, unknown_storage_tag) {
1.1869 + typedef typename E2::size_type size_type;
1.1870 + typedef typename E2::difference_type difference_type;
1.1871 + typedef typename E2::value_type value_type;
1.1872 +
1.1873 + BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
1.1874 + BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
1.1875 + size_type size = e2 ().size ();
1.1876 + for (size_type n = 0; n < size; ++ n) {
1.1877 +#ifndef BOOST_UBLAS_SINGULAR_CHECK
1.1878 + BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
1.1879 +#else
1.1880 + if (e1 () (n, n) == value_type/*zero*/())
1.1881 + singular ().raise ();
1.1882 +#endif
1.1883 + value_type t = e2 () (n) /= e1 () (n, n);
1.1884 + if (t != value_type/*zero*/()) {
1.1885 + typename E1::const_iterator1 it1e1 (e1 ().find1 (1, n + 1, n));
1.1886 + typename E1::const_iterator1 it1e1_end (e1 ().find1 (1, e1 ().size1 (), n));
1.1887 + while (it1e1 != it1e1_end)
1.1888 + e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1;
1.1889 + }
1.1890 + }
1.1891 + }
1.1892 + // Redirectors :-)
1.1893 + template<class E1, class E2>
1.1894 + BOOST_UBLAS_INLINE
1.1895 + void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
1.1896 + lower_tag, column_major_tag) {
1.1897 + typedef typename E1::storage_category storage_category;
1.1898 + inplace_solve (e1, e2,
1.1899 + lower_tag (), column_major_tag (), storage_category ());
1.1900 + }
1.1901 + template<class E1, class E2>
1.1902 + BOOST_UBLAS_INLINE
1.1903 + void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
1.1904 + lower_tag, row_major_tag) {
1.1905 + typedef typename E1::storage_category storage_category;
1.1906 + inplace_solve (e2, trans (e1),
1.1907 + upper_tag (), row_major_tag (), storage_category ());
1.1908 + }
1.1909 + // Dispatcher
1.1910 + template<class E1, class E2>
1.1911 + BOOST_UBLAS_INLINE
1.1912 + void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
1.1913 + lower_tag) {
1.1914 + typedef typename E1::orientation_category orientation_category;
1.1915 + inplace_solve (e1, e2,
1.1916 + lower_tag (), orientation_category ());
1.1917 + }
1.1918 + template<class E1, class E2>
1.1919 + BOOST_UBLAS_INLINE
1.1920 + void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
1.1921 + unit_lower_tag) {
1.1922 + typedef typename E1::orientation_category orientation_category;
1.1923 + inplace_solve (triangular_adaptor<const E1, unit_lower> (e1 ()), e2,
1.1924 + unit_lower_tag (), orientation_category ());
1.1925 + }
1.1926 +
1.1927 + // Dense (proxy) case
1.1928 + template<class E1, class E2>
1.1929 + BOOST_UBLAS_INLINE
1.1930 + void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
1.1931 + upper_tag, column_major_tag, dense_proxy_tag) {
1.1932 + typedef typename E2::size_type size_type;
1.1933 + typedef typename E2::difference_type difference_type;
1.1934 + typedef typename E2::value_type value_type;
1.1935 +
1.1936 + BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
1.1937 + BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
1.1938 + size_type size = e2 ().size ();
1.1939 + for (difference_type n = size - 1; n >= 0; -- n) {
1.1940 +#ifndef BOOST_UBLAS_SINGULAR_CHECK
1.1941 + BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
1.1942 +#else
1.1943 + if (e1 () (n, n) == value_type/*zero*/())
1.1944 + singular ().raise ();
1.1945 +#endif
1.1946 + value_type t = e2 () (n) /= e1 () (n, n);
1.1947 + if (t != value_type/*zero*/()) {
1.1948 + for (difference_type m = n - 1; m >= 0; -- m)
1.1949 + e2 () (m) -= e1 () (m, n) * t;
1.1950 + }
1.1951 + }
1.1952 + }
1.1953 + // Packed (proxy) case
1.1954 + template<class E1, class E2>
1.1955 + BOOST_UBLAS_INLINE
1.1956 + void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
1.1957 + upper_tag, column_major_tag, packed_proxy_tag) {
1.1958 + typedef typename E2::size_type size_type;
1.1959 + typedef typename E2::difference_type difference_type;
1.1960 + typedef typename E2::value_type value_type;
1.1961 +
1.1962 + BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
1.1963 + BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
1.1964 + size_type size = e2 ().size ();
1.1965 + for (difference_type n = size - 1; n >= 0; -- n) {
1.1966 +#ifndef BOOST_UBLAS_SINGULAR_CHECK
1.1967 + BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
1.1968 +#else
1.1969 + if (e1 () (n, n) == value_type/*zero*/())
1.1970 + singular ().raise ();
1.1971 +#endif
1.1972 + value_type t = e2 () (n) /= e1 () (n, n);
1.1973 + if (t != value_type/*zero*/()) {
1.1974 + typename E1::const_reverse_iterator1 it1e1 (e1 ().find1 (1, n, n));
1.1975 + typename E1::const_reverse_iterator1 it1e1_rend (e1 ().find1 (1, 0, n));
1.1976 + difference_type m (it1e1_rend - it1e1);
1.1977 + while (-- m >= 0)
1.1978 + e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1;
1.1979 + }
1.1980 + }
1.1981 + }
1.1982 + // Sparse (proxy) case
1.1983 + template<class E1, class E2>
1.1984 + BOOST_UBLAS_INLINE
1.1985 + void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
1.1986 + upper_tag, column_major_tag, unknown_storage_tag) {
1.1987 + typedef typename E2::size_type size_type;
1.1988 + typedef typename E2::difference_type difference_type;
1.1989 + typedef typename E2::value_type value_type;
1.1990 +
1.1991 + BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
1.1992 + BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ());
1.1993 + size_type size = e2 ().size ();
1.1994 + for (difference_type n = size - 1; n >= 0; -- n) {
1.1995 +#ifndef BOOST_UBLAS_SINGULAR_CHECK
1.1996 + BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
1.1997 +#else
1.1998 + if (e1 () (n, n) == value_type/*zero*/())
1.1999 + singular ().raise ();
1.2000 +#endif
1.2001 + value_type t = e2 () (n) /= e1 () (n, n);
1.2002 + if (t != value_type/*zero*/()) {
1.2003 + typename E1::const_reverse_iterator1 it1e1 (e1 ().find1 (1, n, n));
1.2004 + typename E1::const_reverse_iterator1 it1e1_rend (e1 ().find1 (1, 0, n));
1.2005 + while (it1e1 != it1e1_rend)
1.2006 + e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1;
1.2007 + }
1.2008 + }
1.2009 + }
1.2010 + // Redirectors :-)
1.2011 + template<class E1, class E2>
1.2012 + BOOST_UBLAS_INLINE
1.2013 + void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
1.2014 + upper_tag, column_major_tag) {
1.2015 + typedef typename E1::storage_category storage_category;
1.2016 + inplace_solve (e1, e2,
1.2017 + upper_tag (), column_major_tag (), storage_category ());
1.2018 + }
1.2019 + template<class E1, class E2>
1.2020 + BOOST_UBLAS_INLINE
1.2021 + void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
1.2022 + upper_tag, row_major_tag) {
1.2023 + typedef typename E1::storage_category storage_category;
1.2024 + inplace_solve (e2, trans (e1),
1.2025 + lower_tag (), row_major_tag (), storage_category ());
1.2026 + }
1.2027 + // Dispatcher
1.2028 + template<class E1, class E2>
1.2029 + BOOST_UBLAS_INLINE
1.2030 + void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
1.2031 + upper_tag) {
1.2032 + typedef typename E1::orientation_category orientation_category;
1.2033 + inplace_solve (e1, e2,
1.2034 + upper_tag (), orientation_category ());
1.2035 + }
1.2036 + template<class E1, class E2>
1.2037 + BOOST_UBLAS_INLINE
1.2038 + void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2,
1.2039 + unit_upper_tag) {
1.2040 + typedef typename E1::orientation_category orientation_category;
1.2041 + inplace_solve (triangular_adaptor<const E1, unit_upper> (e1 ()), e2,
1.2042 + unit_upper_tag (), orientation_category ());
1.2043 + }
1.2044 +
1.2045 + template<class E1, class E2, class C>
1.2046 + BOOST_UBLAS_INLINE
1.2047 + typename matrix_vector_solve_traits<E1, E2>::result_type
1.2048 + solve (const matrix_expression<E1> &e1,
1.2049 + const vector_expression<E2> &e2,
1.2050 + C) {
1.2051 + typename matrix_vector_solve_traits<E1, E2>::result_type r (e2);
1.2052 + inplace_solve (e1, r, C ());
1.2053 + return r;
1.2054 + }
1.2055 +
1.2056 + // Dense (proxy) case
1.2057 + template<class E1, class E2>
1.2058 + BOOST_UBLAS_INLINE
1.2059 + void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
1.2060 + lower_tag, row_major_tag, dense_proxy_tag) {
1.2061 + typedef typename E1::size_type size_type;
1.2062 + typedef typename E1::difference_type difference_type;
1.2063 + typedef typename E1::value_type value_type;
1.2064 +
1.2065 + BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());
1.2066 + BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());
1.2067 + size_type size = e1 ().size ();
1.2068 + for (difference_type n = size - 1; n >= 0; -- n) {
1.2069 +#ifndef BOOST_UBLAS_SINGULAR_CHECK
1.2070 + BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());
1.2071 +#else
1.2072 + if (e2 () (n, n) == value_type/*zero*/())
1.2073 + singular ().raise ();
1.2074 +#endif
1.2075 + value_type t = e1 () (n) /= e2 () (n, n);
1.2076 + if (t != value_type/*zero*/()) {
1.2077 + for (difference_type m = n - 1; m >= 0; -- m)
1.2078 + e1 () (m) -= t * e2 () (n, m);
1.2079 + }
1.2080 + }
1.2081 + }
1.2082 + // Packed (proxy) case
1.2083 + template<class E1, class E2>
1.2084 + BOOST_UBLAS_INLINE
1.2085 + void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
1.2086 + lower_tag, row_major_tag, packed_proxy_tag) {
1.2087 + typedef typename E1::size_type size_type;
1.2088 + typedef typename E1::difference_type difference_type;
1.2089 + typedef typename E1::value_type value_type;
1.2090 +
1.2091 + BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());
1.2092 + BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());
1.2093 + size_type size = e1 ().size ();
1.2094 + for (difference_type n = size - 1; n >= 0; -- n) {
1.2095 +#ifndef BOOST_UBLAS_SINGULAR_CHECK
1.2096 + BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());
1.2097 +#else
1.2098 + if (e2 () (n, n) == value_type/*zero*/())
1.2099 + singular ().raise ();
1.2100 +#endif
1.2101 + value_type t = e1 () (n) /= e2 () (n, n);
1.2102 + if (t != value_type/*zero*/()) {
1.2103 + typename E2::const_reverse_iterator2 it2e2 (e2 ().find2 (1, n, n));
1.2104 + typename E2::const_reverse_iterator2 it2e2_rend (e2 ().find2 (1, n, 0));
1.2105 + difference_type m (it2e2_rend - it2e2);
1.2106 + while (-- m >= 0)
1.2107 + e1 () (it2e2.index2 ()) -= *it2e2 * t, ++ it2e2;
1.2108 + }
1.2109 + }
1.2110 + }
1.2111 + // Sparse (proxy) case
1.2112 + template<class E1, class E2>
1.2113 + BOOST_UBLAS_INLINE
1.2114 + void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
1.2115 + lower_tag, row_major_tag, unknown_storage_tag) {
1.2116 + typedef typename E1::size_type size_type;
1.2117 + typedef typename E1::difference_type difference_type;
1.2118 + typedef typename E1::value_type value_type;
1.2119 +
1.2120 + BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());
1.2121 + BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());
1.2122 + size_type size = e1 ().size ();
1.2123 + for (difference_type n = size - 1; n >= 0; -- n) {
1.2124 +#ifndef BOOST_UBLAS_SINGULAR_CHECK
1.2125 + BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());
1.2126 +#else
1.2127 + if (e2 () (n, n) == value_type/*zero*/())
1.2128 + singular ().raise ();
1.2129 +#endif
1.2130 + value_type t = e1 () (n) /= e2 () (n, n);
1.2131 + if (t != value_type/*zero*/()) {
1.2132 + typename E2::const_reverse_iterator2 it2e2 (e2 ().find2 (1, n, n));
1.2133 + typename E2::const_reverse_iterator2 it2e2_rend (e2 ().find2 (1, n, 0));
1.2134 + while (it2e2 != it2e2_rend)
1.2135 + e1 () (it2e2.index2 ()) -= *it2e2 * t, ++ it2e2;
1.2136 + }
1.2137 + }
1.2138 + }
1.2139 + // Redirectors :-)
1.2140 + template<class E1, class E2>
1.2141 + BOOST_UBLAS_INLINE
1.2142 + void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
1.2143 + lower_tag, row_major_tag) {
1.2144 + typedef typename E1::storage_category storage_category;
1.2145 + inplace_solve (e1, e2,
1.2146 + lower_tag (), row_major_tag (), storage_category ());
1.2147 + }
1.2148 + template<class E1, class E2>
1.2149 + BOOST_UBLAS_INLINE
1.2150 + void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
1.2151 + lower_tag, column_major_tag) {
1.2152 + typedef typename E1::storage_category storage_category;
1.2153 + inplace_solve (trans (e2), e1,
1.2154 + upper_tag (), row_major_tag (), storage_category ());
1.2155 + }
1.2156 + // Dispatcher
1.2157 + template<class E1, class E2>
1.2158 + BOOST_UBLAS_INLINE
1.2159 + void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
1.2160 + lower_tag) {
1.2161 + typedef typename E2::orientation_category orientation_category;
1.2162 + inplace_solve (e1, e2,
1.2163 + lower_tag (), orientation_category ());
1.2164 + }
1.2165 + template<class E1, class E2>
1.2166 + BOOST_UBLAS_INLINE
1.2167 + void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
1.2168 + unit_lower_tag) {
1.2169 + typedef typename E2::orientation_category orientation_category;
1.2170 + inplace_solve (e1, triangular_adaptor<const E2, unit_lower> (e2 ()),
1.2171 + unit_lower_tag (), orientation_category ());
1.2172 + }
1.2173 +
1.2174 + // Dense (proxy) case
1.2175 + template<class E1, class E2>
1.2176 + BOOST_UBLAS_INLINE
1.2177 + void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
1.2178 + upper_tag, row_major_tag, dense_proxy_tag) {
1.2179 + typedef typename E1::size_type size_type;
1.2180 + typedef typename E1::difference_type difference_type;
1.2181 + typedef typename E1::value_type value_type;
1.2182 +
1.2183 + BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());
1.2184 + BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());
1.2185 + size_type size = e1 ().size ();
1.2186 + for (size_type n = 0; n < size; ++ n) {
1.2187 +#ifndef BOOST_UBLAS_SINGULAR_CHECK
1.2188 + BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());
1.2189 +#else
1.2190 + if (e2 () (n, n) == value_type/*zero*/())
1.2191 + singular ().raise ();
1.2192 +#endif
1.2193 + value_type t = e1 () (n) /= e2 () (n, n);
1.2194 + if (t != value_type/*zero*/()) {
1.2195 + for (size_type m = n + 1; m < size; ++ m)
1.2196 + e1 () (m) -= t * e2 () (n, m);
1.2197 + }
1.2198 + }
1.2199 + }
1.2200 + // Packed (proxy) case
1.2201 + template<class E1, class E2>
1.2202 + BOOST_UBLAS_INLINE
1.2203 + void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
1.2204 + upper_tag, row_major_tag, packed_proxy_tag) {
1.2205 + typedef typename E1::size_type size_type;
1.2206 + typedef typename E1::difference_type difference_type;
1.2207 + typedef typename E1::value_type value_type;
1.2208 +
1.2209 + BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());
1.2210 + BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());
1.2211 + size_type size = e1 ().size ();
1.2212 + for (size_type n = 0; n < size; ++ n) {
1.2213 +#ifndef BOOST_UBLAS_SINGULAR_CHECK
1.2214 + BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());
1.2215 +#else
1.2216 + if (e2 () (n, n) == value_type/*zero*/())
1.2217 + singular ().raise ();
1.2218 +#endif
1.2219 + value_type t = e1 () (n) /= e2 () (n, n);
1.2220 + if (t != value_type/*zero*/()) {
1.2221 + typename E2::const_iterator2 it2e2 (e2 ().find2 (1, n, n + 1));
1.2222 + typename E2::const_iterator2 it2e2_end (e2 ().find2 (1, n, e2 ().size2 ()));
1.2223 + difference_type m (it2e2_end - it2e2);
1.2224 + while (-- m >= 0)
1.2225 + e1 () (it2e2.index2 ()) -= *it2e2 * t, ++ it2e2;
1.2226 + }
1.2227 + }
1.2228 + }
1.2229 + // Sparse (proxy) case
1.2230 + template<class E1, class E2>
1.2231 + BOOST_UBLAS_INLINE
1.2232 + void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
1.2233 + upper_tag, row_major_tag, unknown_storage_tag) {
1.2234 + typedef typename E1::size_type size_type;
1.2235 + typedef typename E1::difference_type difference_type;
1.2236 + typedef typename E1::value_type value_type;
1.2237 +
1.2238 + BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ());
1.2239 + BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ());
1.2240 + size_type size = e1 ().size ();
1.2241 + for (size_type n = 0; n < size; ++ n) {
1.2242 +#ifndef BOOST_UBLAS_SINGULAR_CHECK
1.2243 + BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ());
1.2244 +#else
1.2245 + if (e2 () (n, n) == value_type/*zero*/())
1.2246 + singular ().raise ();
1.2247 +#endif
1.2248 + value_type t = e1 () (n) /= e2 () (n, n);
1.2249 + if (t != value_type/*zero*/()) {
1.2250 + typename E2::const_iterator2 it2e2 (e2 ().find2 (1, n, n + 1));
1.2251 + typename E2::const_iterator2 it2e2_end (e2 ().find2 (1, n, e2 ().size2 ()));
1.2252 + while (it2e2 != it2e2_end)
1.2253 + e1 () (it2e2.index2 ()) -= *it2e2 * t, ++ it2e2;
1.2254 + }
1.2255 + }
1.2256 + }
1.2257 + // Redirectors :-)
1.2258 + template<class E1, class E2>
1.2259 + BOOST_UBLAS_INLINE
1.2260 + void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
1.2261 + upper_tag, row_major_tag) {
1.2262 + typedef typename E1::storage_category storage_category;
1.2263 + inplace_solve (e1, e2,
1.2264 + upper_tag (), row_major_tag (), storage_category ());
1.2265 + }
1.2266 + template<class E1, class E2>
1.2267 + BOOST_UBLAS_INLINE
1.2268 + void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
1.2269 + upper_tag, column_major_tag) {
1.2270 + typedef typename E1::storage_category storage_category;
1.2271 + inplace_solve (trans (e2), e1,
1.2272 + lower_tag (), row_major_tag (), storage_category ());
1.2273 + }
1.2274 + // Dispatcher
1.2275 + template<class E1, class E2>
1.2276 + BOOST_UBLAS_INLINE
1.2277 + void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
1.2278 + upper_tag) {
1.2279 + typedef typename E2::orientation_category orientation_category;
1.2280 + inplace_solve (e1, e2,
1.2281 + upper_tag (), orientation_category ());
1.2282 + }
1.2283 + template<class E1, class E2>
1.2284 + BOOST_UBLAS_INLINE
1.2285 + void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2,
1.2286 + unit_upper_tag) {
1.2287 + typedef typename E2::orientation_category orientation_category;
1.2288 + inplace_solve (e1, triangular_adaptor<const E2, unit_upper> (e2 ()),
1.2289 + unit_upper_tag (), orientation_category ());
1.2290 + }
1.2291 +
1.2292 + template<class E1, class E2, class C>
1.2293 + BOOST_UBLAS_INLINE
1.2294 + typename matrix_vector_solve_traits<E1, E2>::result_type
1.2295 + solve (const vector_expression<E1> &e1,
1.2296 + const matrix_expression<E2> &e2,
1.2297 + C) {
1.2298 + typename matrix_vector_solve_traits<E1, E2>::result_type r (e1);
1.2299 + inplace_solve (r, e2, C ());
1.2300 + return r;
1.2301 + }
1.2302 +
1.2303 + template<class E1, class E2>
1.2304 + struct matrix_matrix_solve_traits {
1.2305 + typedef typename promote_traits<typename E1::value_type, typename E2::value_type>::promote_type promote_type;
1.2306 + typedef matrix<promote_type> result_type;
1.2307 + };
1.2308 +
1.2309 + // Operations:
1.2310 + // k * n * (n - 1) / 2 + k * n = k * n * (n + 1) / 2 multiplications,
1.2311 + // k * n * (n - 1) / 2 additions
1.2312 +
1.2313 + // Dense (proxy) case
1.2314 + template<class E1, class E2>
1.2315 + BOOST_UBLAS_INLINE
1.2316 + void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
1.2317 + lower_tag, dense_proxy_tag) {
1.2318 + typedef typename E2::size_type size_type;
1.2319 + typedef typename E2::difference_type difference_type;
1.2320 + typedef typename E2::value_type value_type;
1.2321 +
1.2322 + BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
1.2323 + BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());
1.2324 + size_type size1 = e2 ().size1 ();
1.2325 + size_type size2 = e2 ().size2 ();
1.2326 + for (size_type n = 0; n < size1; ++ n) {
1.2327 +#ifndef BOOST_UBLAS_SINGULAR_CHECK
1.2328 + BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
1.2329 +#else
1.2330 + if (e1 () (n, n) == value_type/*zero*/())
1.2331 + singular ().raise ();
1.2332 +#endif
1.2333 + for (size_type l = 0; l < size2; ++ l) {
1.2334 + value_type t = e2 () (n, l) /= e1 () (n, n);
1.2335 + if (t != value_type/*zero*/()) {
1.2336 + for (size_type m = n + 1; m < size1; ++ m)
1.2337 + e2 () (m, l) -= e1 () (m, n) * t;
1.2338 + }
1.2339 + }
1.2340 + }
1.2341 + }
1.2342 + // Packed (proxy) case
1.2343 + template<class E1, class E2>
1.2344 + BOOST_UBLAS_INLINE
1.2345 + void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
1.2346 + lower_tag, packed_proxy_tag) {
1.2347 + typedef typename E2::size_type size_type;
1.2348 + typedef typename E2::difference_type difference_type;
1.2349 + typedef typename E2::value_type value_type;
1.2350 +
1.2351 + BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
1.2352 + BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());
1.2353 + size_type size1 = e2 ().size1 ();
1.2354 + size_type size2 = e2 ().size2 ();
1.2355 + for (size_type n = 0; n < size1; ++ n) {
1.2356 +#ifndef BOOST_UBLAS_SINGULAR_CHECK
1.2357 + BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
1.2358 +#else
1.2359 + if (e1 () (n, n) == value_type/*zero*/())
1.2360 + singular ().raise ();
1.2361 +#endif
1.2362 + for (size_type l = 0; l < size2; ++ l) {
1.2363 + value_type t = e2 () (n, l) /= e1 () (n, n);
1.2364 + if (t != value_type/*zero*/()) {
1.2365 + typename E1::const_iterator1 it1e1 (e1 ().find1 (1, n + 1, n));
1.2366 + typename E1::const_iterator1 it1e1_end (e1 ().find1 (1, e1 ().size1 (), n));
1.2367 + difference_type m (it1e1_end - it1e1);
1.2368 + while (-- m >= 0)
1.2369 + e2 () (it1e1.index1 (), l) -= *it1e1 * t, ++ it1e1;
1.2370 + }
1.2371 + }
1.2372 + }
1.2373 + }
1.2374 + // Sparse (proxy) case
1.2375 + template<class E1, class E2>
1.2376 + BOOST_UBLAS_INLINE
1.2377 + void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
1.2378 + lower_tag, unknown_storage_tag) {
1.2379 + typedef typename E2::size_type size_type;
1.2380 + typedef typename E2::difference_type difference_type;
1.2381 + typedef typename E2::value_type value_type;
1.2382 +
1.2383 + BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
1.2384 + BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());
1.2385 + size_type size1 = e2 ().size1 ();
1.2386 + size_type size2 = e2 ().size2 ();
1.2387 + for (size_type n = 0; n < size1; ++ n) {
1.2388 +#ifndef BOOST_UBLAS_SINGULAR_CHECK
1.2389 + BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
1.2390 +#else
1.2391 + if (e1 () (n, n) == value_type/*zero*/())
1.2392 + singular ().raise ();
1.2393 +#endif
1.2394 + for (size_type l = 0; l < size2; ++ l) {
1.2395 + value_type t = e2 () (n, l) /= e1 () (n, n);
1.2396 + if (t != value_type/*zero*/()) {
1.2397 + typename E1::const_iterator1 it1e1 (e1 ().find1 (1, n + 1, n));
1.2398 + typename E1::const_iterator1 it1e1_end (e1 ().find1 (1, e1 ().size1 (), n));
1.2399 + while (it1e1 != it1e1_end)
1.2400 + e2 () (it1e1.index1 (), l) -= *it1e1 * t, ++ it1e1;
1.2401 + }
1.2402 + }
1.2403 + }
1.2404 + }
1.2405 + // Dispatcher
1.2406 + template<class E1, class E2>
1.2407 + BOOST_UBLAS_INLINE
1.2408 + void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
1.2409 + lower_tag) {
1.2410 + typedef typename E1::storage_category dispatch_category;
1.2411 + inplace_solve (e1, e2,
1.2412 + lower_tag (), dispatch_category ());
1.2413 + }
1.2414 + template<class E1, class E2>
1.2415 + BOOST_UBLAS_INLINE
1.2416 + void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
1.2417 + unit_lower_tag) {
1.2418 + typedef typename E1::storage_category dispatch_category;
1.2419 + inplace_solve (triangular_adaptor<const E1, unit_lower> (e1 ()), e2,
1.2420 + unit_lower_tag (), dispatch_category ());
1.2421 + }
1.2422 +
1.2423 + // Dense (proxy) case
1.2424 + template<class E1, class E2>
1.2425 + BOOST_UBLAS_INLINE
1.2426 + void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
1.2427 + upper_tag, dense_proxy_tag) {
1.2428 + typedef typename E2::size_type size_type;
1.2429 + typedef typename E2::difference_type difference_type;
1.2430 + typedef typename E2::value_type value_type;
1.2431 +
1.2432 + BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
1.2433 + BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());
1.2434 + size_type size1 = e2 ().size1 ();
1.2435 + size_type size2 = e2 ().size2 ();
1.2436 + for (difference_type n = size1 - 1; n >= 0; -- n) {
1.2437 +#ifndef BOOST_UBLAS_SINGULAR_CHECK
1.2438 + BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
1.2439 +#else
1.2440 + if (e1 () (n, n) == value_type/*zero*/())
1.2441 + singular ().raise ();
1.2442 +#endif
1.2443 + for (difference_type l = size2 - 1; l >= 0; -- l) {
1.2444 + value_type t = e2 () (n, l) /= e1 () (n, n);
1.2445 + if (t != value_type/*zero*/()) {
1.2446 + for (difference_type m = n - 1; m >= 0; -- m)
1.2447 + e2 () (m, l) -= e1 () (m, n) * t;
1.2448 + }
1.2449 + }
1.2450 + }
1.2451 + }
1.2452 + // Packed (proxy) case
1.2453 + template<class E1, class E2>
1.2454 + BOOST_UBLAS_INLINE
1.2455 + void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
1.2456 + upper_tag, packed_proxy_tag) {
1.2457 + typedef typename E2::size_type size_type;
1.2458 + typedef typename E2::difference_type difference_type;
1.2459 + typedef typename E2::value_type value_type;
1.2460 +
1.2461 + BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
1.2462 + BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());
1.2463 + size_type size1 = e2 ().size1 ();
1.2464 + size_type size2 = e2 ().size2 ();
1.2465 + for (difference_type n = size1 - 1; n >= 0; -- n) {
1.2466 +#ifndef BOOST_UBLAS_SINGULAR_CHECK
1.2467 + BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
1.2468 +#else
1.2469 + if (e1 () (n, n) == value_type/*zero*/())
1.2470 + singular ().raise ();
1.2471 +#endif
1.2472 + for (difference_type l = size2 - 1; l >= 0; -- l) {
1.2473 + value_type t = e2 () (n, l) /= e1 () (n, n);
1.2474 + if (t != value_type/*zero*/()) {
1.2475 + typename E1::const_reverse_iterator1 it1e1 (e1 ().find1 (1, n, n));
1.2476 + typename E1::const_reverse_iterator1 it1e1_rend (e1 ().find1 (1, 0, n));
1.2477 + difference_type m (it1e1_rend - it1e1);
1.2478 + while (-- m >= 0)
1.2479 + e2 () (it1e1.index1 (), l) -= *it1e1 * t, ++ it1e1;
1.2480 + }
1.2481 + }
1.2482 + }
1.2483 + }
1.2484 + // Sparse (proxy) case
1.2485 + template<class E1, class E2>
1.2486 + BOOST_UBLAS_INLINE
1.2487 + void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
1.2488 + upper_tag, unknown_storage_tag) {
1.2489 + typedef typename E2::size_type size_type;
1.2490 + typedef typename E2::difference_type difference_type;
1.2491 + typedef typename E2::value_type value_type;
1.2492 +
1.2493 + BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
1.2494 + BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ());
1.2495 + size_type size1 = e2 ().size1 ();
1.2496 + size_type size2 = e2 ().size2 ();
1.2497 + for (difference_type n = size1 - 1; n >= 0; -- n) {
1.2498 +#ifndef BOOST_UBLAS_SINGULAR_CHECK
1.2499 + BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ());
1.2500 +#else
1.2501 + if (e1 () (n, n) == value_type/*zero*/())
1.2502 + singular ().raise ();
1.2503 +#endif
1.2504 + for (difference_type l = size2 - 1; l >= 0; -- l) {
1.2505 + value_type t = e2 () (n, l) /= e1 () (n, n);
1.2506 + if (t != value_type/*zero*/()) {
1.2507 + typename E1::const_reverse_iterator1 it1e1 (e1 ().find1 (1, n, n));
1.2508 + typename E1::const_reverse_iterator1 it1e1_rend (e1 ().find1 (1, 0, n));
1.2509 + while (it1e1 != it1e1_rend)
1.2510 + e2 () (it1e1.index1 (), l) -= *it1e1 * t, ++ it1e1;
1.2511 + }
1.2512 + }
1.2513 + }
1.2514 + }
1.2515 + // Dispatcher
1.2516 + template<class E1, class E2>
1.2517 + BOOST_UBLAS_INLINE
1.2518 + void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
1.2519 + upper_tag) {
1.2520 + typedef typename E1::storage_category dispatch_category;
1.2521 + inplace_solve (e1, e2,
1.2522 + upper_tag (), dispatch_category ());
1.2523 + }
1.2524 + template<class E1, class E2>
1.2525 + BOOST_UBLAS_INLINE
1.2526 + void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2,
1.2527 + unit_upper_tag) {
1.2528 + typedef typename E1::storage_category dispatch_category;
1.2529 + inplace_solve (triangular_adaptor<const E1, unit_upper> (e1 ()), e2,
1.2530 + unit_upper_tag (), dispatch_category ());
1.2531 + }
1.2532 +
1.2533 + template<class E1, class E2, class C>
1.2534 + BOOST_UBLAS_INLINE
1.2535 + typename matrix_matrix_solve_traits<E1, E2>::result_type
1.2536 + solve (const matrix_expression<E1> &e1,
1.2537 + const matrix_expression<E2> &e2,
1.2538 + C) {
1.2539 + typename matrix_matrix_solve_traits<E1, E2>::result_type r (e2);
1.2540 + inplace_solve (e1, r, C ());
1.2541 + return r;
1.2542 + }
1.2543 +
1.2544 +}}}
1.2545 +
1.2546 +#endif