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