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