...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
00001 // 00002 // Copyright (c) 2000-2002 00003 // Joerg Walter, Mathias Koch 00004 // 00005 // Distributed under the Boost Software License, Version 1.0. (See 00006 // accompanying file LICENSE_1_0.txt or copy at 00007 // http://www.boost.org/LICENSE_1_0.txt) 00008 // 00009 // The authors gratefully acknowledge the support of 00010 // GeNeSys mbH & Co. KG in producing this work. 00011 // 00012 00013 #ifndef _BOOST_UBLAS_TRIANGULAR_ 00014 #define _BOOST_UBLAS_TRIANGULAR_ 00015 00016 #include <boost/numeric/ublas/matrix.hpp> 00017 #include <boost/numeric/ublas/detail/temporary.hpp> 00018 #include <boost/type_traits/remove_const.hpp> 00019 00020 // Iterators based on ideas of Jeremy Siek 00021 00022 namespace boost { namespace numeric { namespace ublas { 00023 00024 namespace detail { 00025 using namespace boost::numeric::ublas; 00026 00027 // Matrix resizing algorithm 00028 template <class L, class T, class M> 00029 BOOST_UBLAS_INLINE 00030 void matrix_resize_preserve (M& m, M& temporary) { 00031 typedef L layout_type; 00032 typedef T triangular_type; 00033 typedef typename M::size_type size_type; 00034 const size_type msize1 (m.size1 ()); // original size 00035 const size_type msize2 (m.size2 ()); 00036 const size_type size1 (temporary.size1 ()); // new size is specified by temporary 00037 const size_type size2 (temporary.size2 ()); 00038 // Common elements to preserve 00039 const size_type size1_min = (std::min) (size1, msize1); 00040 const size_type size2_min = (std::min) (size2, msize2); 00041 // Order for major and minor sizes 00042 const size_type major_size = layout_type::size_M (size1_min, size2_min); 00043 const size_type minor_size = layout_type::size_m (size1_min, size2_min); 00044 // Indexing copy over major 00045 for (size_type major = 0; major != major_size; ++major) { 00046 for (size_type minor = 0; minor != minor_size; ++minor) { 00047 // find indexes - use invertability of element_ functions 00048 const size_type i1 = layout_type::index_M(major, minor); 00049 const size_type i2 = layout_type::index_m(major, minor); 00050 if ( triangular_type::other(i1,i2) ) { 00051 temporary.data () [triangular_type::element (layout_type (), i1, size1, i2, size2)] = 00052 m.data() [triangular_type::element (layout_type (), i1, msize1, i2, msize2)]; 00053 } 00054 } 00055 } 00056 m.assign_temporary (temporary); 00057 } 00058 } 00059 00077 template<class T, class TRI, class L, class A> 00078 class triangular_matrix: 00079 public matrix_container<triangular_matrix<T, TRI, L, A> > { 00080 00081 typedef T *pointer; 00082 typedef TRI triangular_type; 00083 typedef L layout_type; 00084 typedef triangular_matrix<T, TRI, L, A> self_type; 00085 public: 00086 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS 00087 using matrix_container<self_type>::operator (); 00088 #endif 00089 typedef typename A::size_type size_type; 00090 typedef typename A::difference_type difference_type; 00091 typedef T value_type; 00092 typedef const T &const_reference; 00093 typedef T &reference; 00094 typedef A array_type; 00095 00096 typedef const matrix_reference<const self_type> const_closure_type; 00097 typedef matrix_reference<self_type> closure_type; 00098 typedef vector<T, A> vector_temporary_type; 00099 typedef matrix<T, L, A> matrix_temporary_type; // general sub-matrix 00100 typedef packed_tag storage_category; 00101 typedef typename L::orientation_category orientation_category; 00102 00103 // Construction and destruction 00104 BOOST_UBLAS_INLINE 00105 triangular_matrix (): 00106 matrix_container<self_type> (), 00107 size1_ (0), size2_ (0), data_ (0) {} 00108 BOOST_UBLAS_INLINE 00109 triangular_matrix (size_type size1, size_type size2): 00110 matrix_container<self_type> (), 00111 size1_ (size1), size2_ (size2), data_ (triangular_type::packed_size (layout_type (), size1, size2)) { 00112 } 00113 BOOST_UBLAS_INLINE 00114 triangular_matrix (size_type size1, size_type size2, const array_type &data): 00115 matrix_container<self_type> (), 00116 size1_ (size1), size2_ (size2), data_ (data) {} 00117 BOOST_UBLAS_INLINE 00118 triangular_matrix (const triangular_matrix &m): 00119 matrix_container<self_type> (), 00120 size1_ (m.size1_), size2_ (m.size2_), data_ (m.data_) {} 00121 template<class AE> 00122 BOOST_UBLAS_INLINE 00123 triangular_matrix (const matrix_expression<AE> &ae): 00124 matrix_container<self_type> (), 00125 size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), 00126 data_ (triangular_type::packed_size (layout_type (), size1_, size2_)) { 00127 matrix_assign<scalar_assign> (*this, ae); 00128 } 00129 00130 // Accessors 00131 BOOST_UBLAS_INLINE 00132 size_type size1 () const { 00133 return size1_; 00134 } 00135 BOOST_UBLAS_INLINE 00136 size_type size2 () const { 00137 return size2_; 00138 } 00139 00140 // Storage accessors 00141 BOOST_UBLAS_INLINE 00142 const array_type &data () const { 00143 return data_; 00144 } 00145 BOOST_UBLAS_INLINE 00146 array_type &data () { 00147 return data_; 00148 } 00149 00150 // Resizing 00151 BOOST_UBLAS_INLINE 00152 void resize (size_type size1, size_type size2, bool preserve = true) { 00153 if (preserve) { 00154 self_type temporary (size1, size2); 00155 detail::matrix_resize_preserve<layout_type, triangular_type> (*this, temporary); 00156 } 00157 else { 00158 data ().resize (triangular_type::packed_size (layout_type (), size1, size2)); 00159 size1_ = size1; 00160 size2_ = size2; 00161 } 00162 } 00163 BOOST_UBLAS_INLINE 00164 void resize_packed_preserve (size_type size1, size_type size2) { 00165 size1_ = size1; 00166 size2_ = size2; 00167 data ().resize (triangular_type::packed_size (layout_type (), size1_, size2_), value_type ()); 00168 } 00169 00170 // Element access 00171 BOOST_UBLAS_INLINE 00172 const_reference operator () (size_type i, size_type j) const { 00173 BOOST_UBLAS_CHECK (i < size1_, bad_index ()); 00174 BOOST_UBLAS_CHECK (j < size2_, bad_index ()); 00175 if (triangular_type::other (i, j)) 00176 return data () [triangular_type::element (layout_type (), i, size1_, j, size2_)]; 00177 else if (triangular_type::one (i, j)) 00178 return one_; 00179 else 00180 return zero_; 00181 } 00182 BOOST_UBLAS_INLINE 00183 reference at_element (size_type i, size_type j) { 00184 BOOST_UBLAS_CHECK (i < size1_, bad_index ()); 00185 BOOST_UBLAS_CHECK (j < size2_, bad_index ()); 00186 return data () [triangular_type::element (layout_type (), i, size1_, j, size2_)]; 00187 } 00188 BOOST_UBLAS_INLINE 00189 reference operator () (size_type i, size_type j) { 00190 BOOST_UBLAS_CHECK (i < size1_, bad_index ()); 00191 BOOST_UBLAS_CHECK (j < size2_, bad_index ()); 00192 if (!triangular_type::other (i, j)) { 00193 bad_index ().raise (); 00194 // NEVER reached 00195 } 00196 return data () [triangular_type::element (layout_type (), i, size1_, j, size2_)]; 00197 } 00198 00199 // Element assignment 00200 BOOST_UBLAS_INLINE 00201 reference insert_element (size_type i, size_type j, const_reference t) { 00202 return (operator () (i, j) = t); 00203 } 00204 BOOST_UBLAS_INLINE 00205 void erase_element (size_type i, size_type j) { 00206 operator () (i, j) = value_type/*zero*/(); 00207 } 00208 00209 // Zeroing 00210 BOOST_UBLAS_INLINE 00211 void clear () { 00212 // data ().clear (); 00213 std::fill (data ().begin (), data ().end (), value_type/*zero*/()); 00214 } 00215 00216 // Assignment 00217 BOOST_UBLAS_INLINE 00218 triangular_matrix &operator = (const triangular_matrix &m) { 00219 size1_ = m.size1_; 00220 size2_ = m.size2_; 00221 data () = m.data (); 00222 return *this; 00223 } 00224 BOOST_UBLAS_INLINE 00225 triangular_matrix &assign_temporary (triangular_matrix &m) { 00226 swap (m); 00227 return *this; 00228 } 00229 template<class AE> 00230 BOOST_UBLAS_INLINE 00231 triangular_matrix &operator = (const matrix_expression<AE> &ae) { 00232 self_type temporary (ae); 00233 return assign_temporary (temporary); 00234 } 00235 template<class AE> 00236 BOOST_UBLAS_INLINE 00237 triangular_matrix &assign (const matrix_expression<AE> &ae) { 00238 matrix_assign<scalar_assign> (*this, ae); 00239 return *this; 00240 } 00241 template<class AE> 00242 BOOST_UBLAS_INLINE 00243 triangular_matrix& operator += (const matrix_expression<AE> &ae) { 00244 self_type temporary (*this + ae); 00245 return assign_temporary (temporary); 00246 } 00247 template<class AE> 00248 BOOST_UBLAS_INLINE 00249 triangular_matrix &plus_assign (const matrix_expression<AE> &ae) { 00250 matrix_assign<scalar_plus_assign> (*this, ae); 00251 return *this; 00252 } 00253 template<class AE> 00254 BOOST_UBLAS_INLINE 00255 triangular_matrix& operator -= (const matrix_expression<AE> &ae) { 00256 self_type temporary (*this - ae); 00257 return assign_temporary (temporary); 00258 } 00259 template<class AE> 00260 BOOST_UBLAS_INLINE 00261 triangular_matrix &minus_assign (const matrix_expression<AE> &ae) { 00262 matrix_assign<scalar_minus_assign> (*this, ae); 00263 return *this; 00264 } 00265 template<class AT> 00266 BOOST_UBLAS_INLINE 00267 triangular_matrix& operator *= (const AT &at) { 00268 matrix_assign_scalar<scalar_multiplies_assign> (*this, at); 00269 return *this; 00270 } 00271 template<class AT> 00272 BOOST_UBLAS_INLINE 00273 triangular_matrix& operator /= (const AT &at) { 00274 matrix_assign_scalar<scalar_divides_assign> (*this, at); 00275 return *this; 00276 } 00277 00278 // Swapping 00279 BOOST_UBLAS_INLINE 00280 void swap (triangular_matrix &m) { 00281 if (this != &m) { 00282 // BOOST_UBLAS_CHECK (size2_ == m.size2_, bad_size ()); 00283 std::swap (size1_, m.size1_); 00284 std::swap (size2_, m.size2_); 00285 data ().swap (m.data ()); 00286 } 00287 } 00288 BOOST_UBLAS_INLINE 00289 friend void swap (triangular_matrix &m1, triangular_matrix &m2) { 00290 m1.swap (m2); 00291 } 00292 00293 // Iterator types 00294 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR 00295 typedef indexed_iterator1<self_type, packed_random_access_iterator_tag> iterator1; 00296 typedef indexed_iterator2<self_type, packed_random_access_iterator_tag> iterator2; 00297 typedef indexed_const_iterator1<self_type, packed_random_access_iterator_tag> const_iterator1; 00298 typedef indexed_const_iterator2<self_type, packed_random_access_iterator_tag> const_iterator2; 00299 #else 00300 class const_iterator1; 00301 class iterator1; 00302 class const_iterator2; 00303 class iterator2; 00304 #endif 00305 typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1; 00306 typedef reverse_iterator_base1<iterator1> reverse_iterator1; 00307 typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2; 00308 typedef reverse_iterator_base2<iterator2> reverse_iterator2; 00309 00310 // Element lookup 00311 BOOST_UBLAS_INLINE 00312 const_iterator1 find1 (int rank, size_type i, size_type j) const { 00313 if (rank == 1) 00314 i = triangular_type::restrict1 (i, j, size1_, size2_); 00315 if (rank == 0) 00316 i = triangular_type::global_restrict1 (i, size1_, j, size2_); 00317 return const_iterator1 (*this, i, j); 00318 } 00319 BOOST_UBLAS_INLINE 00320 iterator1 find1 (int rank, size_type i, size_type j) { 00321 if (rank == 1) 00322 i = triangular_type::mutable_restrict1 (i, j, size1_, size2_); 00323 if (rank == 0) 00324 i = triangular_type::global_mutable_restrict1 (i, size1_, j, size2_); 00325 return iterator1 (*this, i, j); 00326 } 00327 BOOST_UBLAS_INLINE 00328 const_iterator2 find2 (int rank, size_type i, size_type j) const { 00329 if (rank == 1) 00330 j = triangular_type::restrict2 (i, j, size1_, size2_); 00331 if (rank == 0) 00332 j = triangular_type::global_restrict2 (i, size1_, j, size2_); 00333 return const_iterator2 (*this, i, j); 00334 } 00335 BOOST_UBLAS_INLINE 00336 iterator2 find2 (int rank, size_type i, size_type j) { 00337 if (rank == 1) 00338 j = triangular_type::mutable_restrict2 (i, j, size1_, size2_); 00339 if (rank == 0) 00340 j = triangular_type::global_mutable_restrict2 (i, size1_, j, size2_); 00341 return iterator2 (*this, i, j); 00342 } 00343 00344 // Iterators simply are indices. 00345 00346 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR 00347 class const_iterator1: 00348 public container_const_reference<triangular_matrix>, 00349 public random_access_iterator_base<packed_random_access_iterator_tag, 00350 const_iterator1, value_type> { 00351 public: 00352 typedef typename triangular_matrix::value_type value_type; 00353 typedef typename triangular_matrix::difference_type difference_type; 00354 typedef typename triangular_matrix::const_reference reference; 00355 typedef const typename triangular_matrix::pointer pointer; 00356 00357 typedef const_iterator2 dual_iterator_type; 00358 typedef const_reverse_iterator2 dual_reverse_iterator_type; 00359 00360 // Construction and destruction 00361 BOOST_UBLAS_INLINE 00362 const_iterator1 (): 00363 container_const_reference<self_type> (), it1_ (), it2_ () {} 00364 BOOST_UBLAS_INLINE 00365 const_iterator1 (const self_type &m, size_type it1, size_type it2): 00366 container_const_reference<self_type> (m), it1_ (it1), it2_ (it2) {} 00367 BOOST_UBLAS_INLINE 00368 const_iterator1 (const iterator1 &it): 00369 container_const_reference<self_type> (it ()), it1_ (it.it1_), it2_ (it.it2_) {} 00370 00371 // Arithmetic 00372 BOOST_UBLAS_INLINE 00373 const_iterator1 &operator ++ () { 00374 ++ it1_; 00375 return *this; 00376 } 00377 BOOST_UBLAS_INLINE 00378 const_iterator1 &operator -- () { 00379 -- it1_; 00380 return *this; 00381 } 00382 BOOST_UBLAS_INLINE 00383 const_iterator1 &operator += (difference_type n) { 00384 it1_ += n; 00385 return *this; 00386 } 00387 BOOST_UBLAS_INLINE 00388 const_iterator1 &operator -= (difference_type n) { 00389 it1_ -= n; 00390 return *this; 00391 } 00392 BOOST_UBLAS_INLINE 00393 difference_type operator - (const const_iterator1 &it) const { 00394 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 00395 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ()); 00396 return it1_ - it.it1_; 00397 } 00398 00399 // Dereference 00400 BOOST_UBLAS_INLINE 00401 const_reference operator * () const { 00402 return (*this) () (it1_, it2_); 00403 } 00404 BOOST_UBLAS_INLINE 00405 const_reference operator [] (difference_type n) const { 00406 return *(*this + n); 00407 } 00408 00409 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 00410 BOOST_UBLAS_INLINE 00411 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 00412 typename self_type:: 00413 #endif 00414 const_iterator2 begin () const { 00415 return (*this) ().find2 (1, it1_, 0); 00416 } 00417 BOOST_UBLAS_INLINE 00418 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 00419 typename self_type:: 00420 #endif 00421 const_iterator2 end () const { 00422 return (*this) ().find2 (1, it1_, (*this) ().size2 ()); 00423 } 00424 BOOST_UBLAS_INLINE 00425 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 00426 typename self_type:: 00427 #endif 00428 const_reverse_iterator2 rbegin () const { 00429 return const_reverse_iterator2 (end ()); 00430 } 00431 BOOST_UBLAS_INLINE 00432 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 00433 typename self_type:: 00434 #endif 00435 const_reverse_iterator2 rend () const { 00436 return const_reverse_iterator2 (begin ()); 00437 } 00438 #endif 00439 00440 // Indices 00441 BOOST_UBLAS_INLINE 00442 size_type index1 () const { 00443 return it1_; 00444 } 00445 BOOST_UBLAS_INLINE 00446 size_type index2 () const { 00447 return it2_; 00448 } 00449 00450 // Assignment 00451 BOOST_UBLAS_INLINE 00452 const_iterator1 &operator = (const const_iterator1 &it) { 00453 container_const_reference<self_type>::assign (&it ()); 00454 it1_ = it.it1_; 00455 it2_ = it.it2_; 00456 return *this; 00457 } 00458 00459 // Comparison 00460 BOOST_UBLAS_INLINE 00461 bool operator == (const const_iterator1 &it) const { 00462 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 00463 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ()); 00464 return it1_ == it.it1_; 00465 } 00466 BOOST_UBLAS_INLINE 00467 bool operator < (const const_iterator1 &it) const { 00468 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 00469 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ()); 00470 return it1_ < it.it1_; 00471 } 00472 00473 private: 00474 size_type it1_; 00475 size_type it2_; 00476 }; 00477 #endif 00478 00479 BOOST_UBLAS_INLINE 00480 const_iterator1 begin1 () const { 00481 return find1 (0, 0, 0); 00482 } 00483 BOOST_UBLAS_INLINE 00484 const_iterator1 end1 () const { 00485 return find1 (0, size1_, 0); 00486 } 00487 00488 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR 00489 class iterator1: 00490 public container_reference<triangular_matrix>, 00491 public random_access_iterator_base<packed_random_access_iterator_tag, 00492 iterator1, value_type> { 00493 public: 00494 typedef typename triangular_matrix::value_type value_type; 00495 typedef typename triangular_matrix::difference_type difference_type; 00496 typedef typename triangular_matrix::reference reference; 00497 typedef typename triangular_matrix::pointer pointer; 00498 00499 typedef iterator2 dual_iterator_type; 00500 typedef reverse_iterator2 dual_reverse_iterator_type; 00501 00502 // Construction and destruction 00503 BOOST_UBLAS_INLINE 00504 iterator1 (): 00505 container_reference<self_type> (), it1_ (), it2_ () {} 00506 BOOST_UBLAS_INLINE 00507 iterator1 (self_type &m, size_type it1, size_type it2): 00508 container_reference<self_type> (m), it1_ (it1), it2_ (it2) {} 00509 00510 // Arithmetic 00511 BOOST_UBLAS_INLINE 00512 iterator1 &operator ++ () { 00513 ++ it1_; 00514 return *this; 00515 } 00516 BOOST_UBLAS_INLINE 00517 iterator1 &operator -- () { 00518 -- it1_; 00519 return *this; 00520 } 00521 BOOST_UBLAS_INLINE 00522 iterator1 &operator += (difference_type n) { 00523 it1_ += n; 00524 return *this; 00525 } 00526 BOOST_UBLAS_INLINE 00527 iterator1 &operator -= (difference_type n) { 00528 it1_ -= n; 00529 return *this; 00530 } 00531 BOOST_UBLAS_INLINE 00532 difference_type operator - (const iterator1 &it) const { 00533 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 00534 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ()); 00535 return it1_ - it.it1_; 00536 } 00537 00538 // Dereference 00539 BOOST_UBLAS_INLINE 00540 reference operator * () const { 00541 return (*this) () (it1_, it2_); 00542 } 00543 BOOST_UBLAS_INLINE 00544 reference operator [] (difference_type n) const { 00545 return *(*this + n); 00546 } 00547 00548 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 00549 BOOST_UBLAS_INLINE 00550 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 00551 typename self_type:: 00552 #endif 00553 iterator2 begin () const { 00554 return (*this) ().find2 (1, it1_, 0); 00555 } 00556 BOOST_UBLAS_INLINE 00557 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 00558 typename self_type:: 00559 #endif 00560 iterator2 end () const { 00561 return (*this) ().find2 (1, it1_, (*this) ().size2 ()); 00562 } 00563 BOOST_UBLAS_INLINE 00564 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 00565 typename self_type:: 00566 #endif 00567 reverse_iterator2 rbegin () const { 00568 return reverse_iterator2 (end ()); 00569 } 00570 BOOST_UBLAS_INLINE 00571 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 00572 typename self_type:: 00573 #endif 00574 reverse_iterator2 rend () const { 00575 return reverse_iterator2 (begin ()); 00576 } 00577 #endif 00578 00579 // Indices 00580 BOOST_UBLAS_INLINE 00581 size_type index1 () const { 00582 return it1_; 00583 } 00584 BOOST_UBLAS_INLINE 00585 size_type index2 () const { 00586 return it2_; 00587 } 00588 00589 // Assignment 00590 BOOST_UBLAS_INLINE 00591 iterator1 &operator = (const iterator1 &it) { 00592 container_reference<self_type>::assign (&it ()); 00593 it1_ = it.it1_; 00594 it2_ = it.it2_; 00595 return *this; 00596 } 00597 00598 // Comparison 00599 BOOST_UBLAS_INLINE 00600 bool operator == (const iterator1 &it) const { 00601 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 00602 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ()); 00603 return it1_ == it.it1_; 00604 } 00605 BOOST_UBLAS_INLINE 00606 bool operator < (const iterator1 &it) const { 00607 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 00608 BOOST_UBLAS_CHECK (it2_ == it.it2_, external_logic ()); 00609 return it1_ < it.it1_; 00610 } 00611 00612 private: 00613 size_type it1_; 00614 size_type it2_; 00615 00616 friend class const_iterator1; 00617 }; 00618 #endif 00619 00620 BOOST_UBLAS_INLINE 00621 iterator1 begin1 () { 00622 return find1 (0, 0, 0); 00623 } 00624 BOOST_UBLAS_INLINE 00625 iterator1 end1 () { 00626 return find1 (0, size1_, 0); 00627 } 00628 00629 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR 00630 class const_iterator2: 00631 public container_const_reference<triangular_matrix>, 00632 public random_access_iterator_base<packed_random_access_iterator_tag, 00633 const_iterator2, value_type> { 00634 public: 00635 typedef typename triangular_matrix::value_type value_type; 00636 typedef typename triangular_matrix::difference_type difference_type; 00637 typedef typename triangular_matrix::const_reference reference; 00638 typedef const typename triangular_matrix::pointer pointer; 00639 00640 typedef const_iterator1 dual_iterator_type; 00641 typedef const_reverse_iterator1 dual_reverse_iterator_type; 00642 00643 // Construction and destruction 00644 BOOST_UBLAS_INLINE 00645 const_iterator2 (): 00646 container_const_reference<self_type> (), it1_ (), it2_ () {} 00647 BOOST_UBLAS_INLINE 00648 const_iterator2 (const self_type &m, size_type it1, size_type it2): 00649 container_const_reference<self_type> (m), it1_ (it1), it2_ (it2) {} 00650 BOOST_UBLAS_INLINE 00651 const_iterator2 (const iterator2 &it): 00652 container_const_reference<self_type> (it ()), it1_ (it.it1_), it2_ (it.it2_) {} 00653 00654 // Arithmetic 00655 BOOST_UBLAS_INLINE 00656 const_iterator2 &operator ++ () { 00657 ++ it2_; 00658 return *this; 00659 } 00660 BOOST_UBLAS_INLINE 00661 const_iterator2 &operator -- () { 00662 -- it2_; 00663 return *this; 00664 } 00665 BOOST_UBLAS_INLINE 00666 const_iterator2 &operator += (difference_type n) { 00667 it2_ += n; 00668 return *this; 00669 } 00670 BOOST_UBLAS_INLINE 00671 const_iterator2 &operator -= (difference_type n) { 00672 it2_ -= n; 00673 return *this; 00674 } 00675 BOOST_UBLAS_INLINE 00676 difference_type operator - (const const_iterator2 &it) const { 00677 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 00678 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ()); 00679 return it2_ - it.it2_; 00680 } 00681 00682 // Dereference 00683 BOOST_UBLAS_INLINE 00684 const_reference operator * () const { 00685 return (*this) () (it1_, it2_); 00686 } 00687 BOOST_UBLAS_INLINE 00688 const_reference operator [] (difference_type n) const { 00689 return *(*this + n); 00690 } 00691 00692 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 00693 BOOST_UBLAS_INLINE 00694 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 00695 typename self_type:: 00696 #endif 00697 const_iterator1 begin () const { 00698 return (*this) ().find1 (1, 0, it2_); 00699 } 00700 BOOST_UBLAS_INLINE 00701 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 00702 typename self_type:: 00703 #endif 00704 const_iterator1 end () const { 00705 return (*this) ().find1 (1, (*this) ().size1 (), it2_); 00706 } 00707 BOOST_UBLAS_INLINE 00708 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 00709 typename self_type:: 00710 #endif 00711 const_reverse_iterator1 rbegin () const { 00712 return const_reverse_iterator1 (end ()); 00713 } 00714 BOOST_UBLAS_INLINE 00715 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 00716 typename self_type:: 00717 #endif 00718 const_reverse_iterator1 rend () const { 00719 return const_reverse_iterator1 (begin ()); 00720 } 00721 #endif 00722 00723 // Indices 00724 BOOST_UBLAS_INLINE 00725 size_type index1 () const { 00726 return it1_; 00727 } 00728 BOOST_UBLAS_INLINE 00729 size_type index2 () const { 00730 return it2_; 00731 } 00732 00733 // Assignment 00734 BOOST_UBLAS_INLINE 00735 const_iterator2 &operator = (const const_iterator2 &it) { 00736 container_const_reference<self_type>::assign (&it ()); 00737 it1_ = it.it1_; 00738 it2_ = it.it2_; 00739 return *this; 00740 } 00741 00742 // Comparison 00743 BOOST_UBLAS_INLINE 00744 bool operator == (const const_iterator2 &it) const { 00745 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 00746 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ()); 00747 return it2_ == it.it2_; 00748 } 00749 BOOST_UBLAS_INLINE 00750 bool operator < (const const_iterator2 &it) const { 00751 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 00752 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ()); 00753 return it2_ < it.it2_; 00754 } 00755 00756 private: 00757 size_type it1_; 00758 size_type it2_; 00759 }; 00760 #endif 00761 00762 BOOST_UBLAS_INLINE 00763 const_iterator2 begin2 () const { 00764 return find2 (0, 0, 0); 00765 } 00766 BOOST_UBLAS_INLINE 00767 const_iterator2 end2 () const { 00768 return find2 (0, 0, size2_); 00769 } 00770 00771 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR 00772 class iterator2: 00773 public container_reference<triangular_matrix>, 00774 public random_access_iterator_base<packed_random_access_iterator_tag, 00775 iterator2, value_type> { 00776 public: 00777 typedef typename triangular_matrix::value_type value_type; 00778 typedef typename triangular_matrix::difference_type difference_type; 00779 typedef typename triangular_matrix::reference reference; 00780 typedef typename triangular_matrix::pointer pointer; 00781 00782 typedef iterator1 dual_iterator_type; 00783 typedef reverse_iterator1 dual_reverse_iterator_type; 00784 00785 // Construction and destruction 00786 BOOST_UBLAS_INLINE 00787 iterator2 (): 00788 container_reference<self_type> (), it1_ (), it2_ () {} 00789 BOOST_UBLAS_INLINE 00790 iterator2 (self_type &m, size_type it1, size_type it2): 00791 container_reference<self_type> (m), it1_ (it1), it2_ (it2) {} 00792 00793 // Arithmetic 00794 BOOST_UBLAS_INLINE 00795 iterator2 &operator ++ () { 00796 ++ it2_; 00797 return *this; 00798 } 00799 BOOST_UBLAS_INLINE 00800 iterator2 &operator -- () { 00801 -- it2_; 00802 return *this; 00803 } 00804 BOOST_UBLAS_INLINE 00805 iterator2 &operator += (difference_type n) { 00806 it2_ += n; 00807 return *this; 00808 } 00809 BOOST_UBLAS_INLINE 00810 iterator2 &operator -= (difference_type n) { 00811 it2_ -= n; 00812 return *this; 00813 } 00814 BOOST_UBLAS_INLINE 00815 difference_type operator - (const iterator2 &it) const { 00816 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 00817 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ()); 00818 return it2_ - it.it2_; 00819 } 00820 00821 // Dereference 00822 BOOST_UBLAS_INLINE 00823 reference operator * () const { 00824 return (*this) () (it1_, it2_); 00825 } 00826 BOOST_UBLAS_INLINE 00827 reference operator [] (difference_type n) const { 00828 return *(*this + n); 00829 } 00830 00831 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 00832 BOOST_UBLAS_INLINE 00833 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 00834 typename self_type:: 00835 #endif 00836 iterator1 begin () const { 00837 return (*this) ().find1 (1, 0, it2_); 00838 } 00839 BOOST_UBLAS_INLINE 00840 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 00841 typename self_type:: 00842 #endif 00843 iterator1 end () const { 00844 return (*this) ().find1 (1, (*this) ().size1 (), it2_); 00845 } 00846 BOOST_UBLAS_INLINE 00847 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 00848 typename self_type:: 00849 #endif 00850 reverse_iterator1 rbegin () const { 00851 return reverse_iterator1 (end ()); 00852 } 00853 BOOST_UBLAS_INLINE 00854 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 00855 typename self_type:: 00856 #endif 00857 reverse_iterator1 rend () const { 00858 return reverse_iterator1 (begin ()); 00859 } 00860 #endif 00861 00862 // Indices 00863 BOOST_UBLAS_INLINE 00864 size_type index1 () const { 00865 return it1_; 00866 } 00867 BOOST_UBLAS_INLINE 00868 size_type index2 () const { 00869 return it2_; 00870 } 00871 00872 // Assignment 00873 BOOST_UBLAS_INLINE 00874 iterator2 &operator = (const iterator2 &it) { 00875 container_reference<self_type>::assign (&it ()); 00876 it1_ = it.it1_; 00877 it2_ = it.it2_; 00878 return *this; 00879 } 00880 00881 // Comparison 00882 BOOST_UBLAS_INLINE 00883 bool operator == (const iterator2 &it) const { 00884 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 00885 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ()); 00886 return it2_ == it.it2_; 00887 } 00888 BOOST_UBLAS_INLINE 00889 bool operator < (const iterator2 &it) const { 00890 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 00891 BOOST_UBLAS_CHECK (it1_ == it.it1_, external_logic ()); 00892 return it2_ < it.it2_; 00893 } 00894 00895 private: 00896 size_type it1_; 00897 size_type it2_; 00898 00899 friend class const_iterator2; 00900 }; 00901 #endif 00902 00903 BOOST_UBLAS_INLINE 00904 iterator2 begin2 () { 00905 return find2 (0, 0, 0); 00906 } 00907 BOOST_UBLAS_INLINE 00908 iterator2 end2 () { 00909 return find2 (0, 0, size2_); 00910 } 00911 00912 // Reverse iterators 00913 00914 BOOST_UBLAS_INLINE 00915 const_reverse_iterator1 rbegin1 () const { 00916 return const_reverse_iterator1 (end1 ()); 00917 } 00918 BOOST_UBLAS_INLINE 00919 const_reverse_iterator1 rend1 () const { 00920 return const_reverse_iterator1 (begin1 ()); 00921 } 00922 00923 BOOST_UBLAS_INLINE 00924 reverse_iterator1 rbegin1 () { 00925 return reverse_iterator1 (end1 ()); 00926 } 00927 BOOST_UBLAS_INLINE 00928 reverse_iterator1 rend1 () { 00929 return reverse_iterator1 (begin1 ()); 00930 } 00931 00932 BOOST_UBLAS_INLINE 00933 const_reverse_iterator2 rbegin2 () const { 00934 return const_reverse_iterator2 (end2 ()); 00935 } 00936 BOOST_UBLAS_INLINE 00937 const_reverse_iterator2 rend2 () const { 00938 return const_reverse_iterator2 (begin2 ()); 00939 } 00940 00941 BOOST_UBLAS_INLINE 00942 reverse_iterator2 rbegin2 () { 00943 return reverse_iterator2 (end2 ()); 00944 } 00945 BOOST_UBLAS_INLINE 00946 reverse_iterator2 rend2 () { 00947 return reverse_iterator2 (begin2 ()); 00948 } 00949 00950 private: 00951 size_type size1_; 00952 size_type size2_; 00953 array_type data_; 00954 static const value_type zero_; 00955 static const value_type one_; 00956 }; 00957 00958 template<class T, class TRI, class L, class A> 00959 const typename triangular_matrix<T, TRI, L, A>::value_type triangular_matrix<T, TRI, L, A>::zero_ = value_type/*zero*/(); 00960 template<class T, class TRI, class L, class A> 00961 const typename triangular_matrix<T, TRI, L, A>::value_type triangular_matrix<T, TRI, L, A>::one_ (1); 00962 00963 00964 // Triangular matrix adaptor class 00965 template<class M, class TRI> 00966 class triangular_adaptor: 00967 public matrix_expression<triangular_adaptor<M, TRI> > { 00968 00969 typedef triangular_adaptor<M, TRI> self_type; 00970 00971 public: 00972 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS 00973 using matrix_expression<self_type>::operator (); 00974 #endif 00975 typedef const M const_matrix_type; 00976 typedef M matrix_type; 00977 typedef TRI triangular_type; 00978 typedef typename M::size_type size_type; 00979 typedef typename M::difference_type difference_type; 00980 typedef typename M::value_type value_type; 00981 typedef typename M::const_reference const_reference; 00982 typedef typename boost::mpl::if_<boost::is_const<M>, 00983 typename M::const_reference, 00984 typename M::reference>::type reference; 00985 typedef typename boost::mpl::if_<boost::is_const<M>, 00986 typename M::const_closure_type, 00987 typename M::closure_type>::type matrix_closure_type; 00988 typedef const self_type const_closure_type; 00989 typedef self_type closure_type; 00990 // Replaced by _temporary_traits to avoid type requirements on M 00991 //typedef typename M::vector_temporary_type vector_temporary_type; 00992 //typedef typename M::matrix_temporary_type matrix_temporary_type; 00993 typedef typename storage_restrict_traits<typename M::storage_category, 00994 packed_proxy_tag>::storage_category storage_category; 00995 typedef typename M::orientation_category orientation_category; 00996 00997 // Construction and destruction 00998 BOOST_UBLAS_INLINE 00999 triangular_adaptor (matrix_type &data): 01000 matrix_expression<self_type> (), 01001 data_ (data) {} 01002 BOOST_UBLAS_INLINE 01003 triangular_adaptor (const triangular_adaptor &m): 01004 matrix_expression<self_type> (), 01005 data_ (m.data_) {} 01006 01007 // Accessors 01008 BOOST_UBLAS_INLINE 01009 size_type size1 () const { 01010 return data_.size1 (); 01011 } 01012 BOOST_UBLAS_INLINE 01013 size_type size2 () const { 01014 return data_.size2 (); 01015 } 01016 01017 // Storage accessors 01018 BOOST_UBLAS_INLINE 01019 const matrix_closure_type &data () const { 01020 return data_; 01021 } 01022 BOOST_UBLAS_INLINE 01023 matrix_closure_type &data () { 01024 return data_; 01025 } 01026 01027 // Element access 01028 #ifndef BOOST_UBLAS_PROXY_CONST_MEMBER 01029 BOOST_UBLAS_INLINE 01030 const_reference operator () (size_type i, size_type j) const { 01031 BOOST_UBLAS_CHECK (i < size1 (), bad_index ()); 01032 BOOST_UBLAS_CHECK (j < size2 (), bad_index ()); 01033 if (triangular_type::other (i, j)) 01034 return data () (i, j); 01035 else if (triangular_type::one (i, j)) 01036 return one_; 01037 else 01038 return zero_; 01039 } 01040 BOOST_UBLAS_INLINE 01041 reference operator () (size_type i, size_type j) { 01042 BOOST_UBLAS_CHECK (i < size1 (), bad_index ()); 01043 BOOST_UBLAS_CHECK (j < size2 (), bad_index ()); 01044 if (!triangular_type::other (i, j)) { 01045 bad_index ().raise (); 01046 // NEVER reached 01047 } 01048 return data () (i, j); 01049 } 01050 #else 01051 BOOST_UBLAS_INLINE 01052 reference operator () (size_type i, size_type j) const { 01053 BOOST_UBLAS_CHECK (i < size1 (), bad_index ()); 01054 BOOST_UBLAS_CHECK (j < size2 (), bad_index ()); 01055 if (!triangular_type::other (i, j)) { 01056 bad_index ().raise (); 01057 // NEVER reached 01058 } 01059 return data () (i, j); 01060 } 01061 #endif 01062 01063 // Assignment 01064 BOOST_UBLAS_INLINE 01065 triangular_adaptor &operator = (const triangular_adaptor &m) { 01066 matrix_assign<scalar_assign> (*this, m); 01067 return *this; 01068 } 01069 BOOST_UBLAS_INLINE 01070 triangular_adaptor &assign_temporary (triangular_adaptor &m) { 01071 *this = m; 01072 return *this; 01073 } 01074 template<class AE> 01075 BOOST_UBLAS_INLINE 01076 triangular_adaptor &operator = (const matrix_expression<AE> &ae) { 01077 matrix_assign<scalar_assign> (*this, matrix<value_type> (ae)); 01078 return *this; 01079 } 01080 template<class AE> 01081 BOOST_UBLAS_INLINE 01082 triangular_adaptor &assign (const matrix_expression<AE> &ae) { 01083 matrix_assign<scalar_assign> (*this, ae); 01084 return *this; 01085 } 01086 template<class AE> 01087 BOOST_UBLAS_INLINE 01088 triangular_adaptor& operator += (const matrix_expression<AE> &ae) { 01089 matrix_assign<scalar_assign> (*this, matrix<value_type> (*this + ae)); 01090 return *this; 01091 } 01092 template<class AE> 01093 BOOST_UBLAS_INLINE 01094 triangular_adaptor &plus_assign (const matrix_expression<AE> &ae) { 01095 matrix_assign<scalar_plus_assign> (*this, ae); 01096 return *this; 01097 } 01098 template<class AE> 01099 BOOST_UBLAS_INLINE 01100 triangular_adaptor& operator -= (const matrix_expression<AE> &ae) { 01101 matrix_assign<scalar_assign> (*this, matrix<value_type> (*this - ae)); 01102 return *this; 01103 } 01104 template<class AE> 01105 BOOST_UBLAS_INLINE 01106 triangular_adaptor &minus_assign (const matrix_expression<AE> &ae) { 01107 matrix_assign<scalar_minus_assign> (*this, ae); 01108 return *this; 01109 } 01110 template<class AT> 01111 BOOST_UBLAS_INLINE 01112 triangular_adaptor& operator *= (const AT &at) { 01113 matrix_assign_scalar<scalar_multiplies_assign> (*this, at); 01114 return *this; 01115 } 01116 template<class AT> 01117 BOOST_UBLAS_INLINE 01118 triangular_adaptor& operator /= (const AT &at) { 01119 matrix_assign_scalar<scalar_divides_assign> (*this, at); 01120 return *this; 01121 } 01122 01123 // Closure comparison 01124 BOOST_UBLAS_INLINE 01125 bool same_closure (const triangular_adaptor &ta) const { 01126 return (*this).data ().same_closure (ta.data ()); 01127 } 01128 01129 // Swapping 01130 BOOST_UBLAS_INLINE 01131 void swap (triangular_adaptor &m) { 01132 if (this != &m) 01133 matrix_swap<scalar_swap> (*this, m); 01134 } 01135 BOOST_UBLAS_INLINE 01136 friend void swap (triangular_adaptor &m1, triangular_adaptor &m2) { 01137 m1.swap (m2); 01138 } 01139 01140 // Iterator types 01141 private: 01142 typedef typename M::const_iterator1 const_subiterator1_type; 01143 typedef typename boost::mpl::if_<boost::is_const<M>, 01144 typename M::const_iterator1, 01145 typename M::iterator1>::type subiterator1_type; 01146 typedef typename M::const_iterator2 const_subiterator2_type; 01147 typedef typename boost::mpl::if_<boost::is_const<M>, 01148 typename M::const_iterator2, 01149 typename M::iterator2>::type subiterator2_type; 01150 01151 public: 01152 #ifdef BOOST_UBLAS_USE_INDEXED_ITERATOR 01153 typedef indexed_iterator1<self_type, packed_random_access_iterator_tag> iterator1; 01154 typedef indexed_iterator2<self_type, packed_random_access_iterator_tag> iterator2; 01155 typedef indexed_const_iterator1<self_type, packed_random_access_iterator_tag> const_iterator1; 01156 typedef indexed_const_iterator2<self_type, packed_random_access_iterator_tag> const_iterator2; 01157 #else 01158 class const_iterator1; 01159 class iterator1; 01160 class const_iterator2; 01161 class iterator2; 01162 #endif 01163 typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1; 01164 typedef reverse_iterator_base1<iterator1> reverse_iterator1; 01165 typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2; 01166 typedef reverse_iterator_base2<iterator2> reverse_iterator2; 01167 01168 // Element lookup 01169 BOOST_UBLAS_INLINE 01170 const_iterator1 find1 (int rank, size_type i, size_type j) const { 01171 if (rank == 1) 01172 i = triangular_type::restrict1 (i, j, size1(), size2()); 01173 if (rank == 0) 01174 i = triangular_type::global_restrict1 (i, size1(), j, size2()); 01175 return const_iterator1 (*this, data ().find1 (rank, i, j)); 01176 } 01177 BOOST_UBLAS_INLINE 01178 iterator1 find1 (int rank, size_type i, size_type j) { 01179 if (rank == 1) 01180 i = triangular_type::mutable_restrict1 (i, j, size1(), size2()); 01181 if (rank == 0) 01182 i = triangular_type::global_mutable_restrict1 (i, size1(), j, size2()); 01183 return iterator1 (*this, data ().find1 (rank, i, j)); 01184 } 01185 BOOST_UBLAS_INLINE 01186 const_iterator2 find2 (int rank, size_type i, size_type j) const { 01187 if (rank == 1) 01188 j = triangular_type::restrict2 (i, j, size1(), size2()); 01189 if (rank == 0) 01190 j = triangular_type::global_restrict2 (i, size1(), j, size2()); 01191 return const_iterator2 (*this, data ().find2 (rank, i, j)); 01192 } 01193 BOOST_UBLAS_INLINE 01194 iterator2 find2 (int rank, size_type i, size_type j) { 01195 if (rank == 1) 01196 j = triangular_type::mutable_restrict2 (i, j, size1(), size2()); 01197 if (rank == 0) 01198 j = triangular_type::global_mutable_restrict2 (i, size1(), j, size2()); 01199 return iterator2 (*this, data ().find2 (rank, i, j)); 01200 } 01201 01202 // Iterators simply are indices. 01203 01204 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR 01205 class const_iterator1: 01206 public container_const_reference<triangular_adaptor>, 01207 public random_access_iterator_base<typename iterator_restrict_traits< 01208 typename const_subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category, 01209 const_iterator1, value_type> { 01210 public: 01211 typedef typename const_subiterator1_type::value_type value_type; 01212 typedef typename const_subiterator1_type::difference_type difference_type; 01213 typedef typename const_subiterator1_type::reference reference; 01214 typedef typename const_subiterator1_type::pointer pointer; 01215 01216 typedef const_iterator2 dual_iterator_type; 01217 typedef const_reverse_iterator2 dual_reverse_iterator_type; 01218 01219 // Construction and destruction 01220 BOOST_UBLAS_INLINE 01221 const_iterator1 (): 01222 container_const_reference<self_type> (), it1_ () {} 01223 BOOST_UBLAS_INLINE 01224 const_iterator1 (const self_type &m, const const_subiterator1_type &it1): 01225 container_const_reference<self_type> (m), it1_ (it1) {} 01226 BOOST_UBLAS_INLINE 01227 const_iterator1 (const iterator1 &it): 01228 container_const_reference<self_type> (it ()), it1_ (it.it1_) {} 01229 01230 // Arithmetic 01231 BOOST_UBLAS_INLINE 01232 const_iterator1 &operator ++ () { 01233 ++ it1_; 01234 return *this; 01235 } 01236 BOOST_UBLAS_INLINE 01237 const_iterator1 &operator -- () { 01238 -- it1_; 01239 return *this; 01240 } 01241 BOOST_UBLAS_INLINE 01242 const_iterator1 &operator += (difference_type n) { 01243 it1_ += n; 01244 return *this; 01245 } 01246 BOOST_UBLAS_INLINE 01247 const_iterator1 &operator -= (difference_type n) { 01248 it1_ -= n; 01249 return *this; 01250 } 01251 BOOST_UBLAS_INLINE 01252 difference_type operator - (const const_iterator1 &it) const { 01253 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 01254 return it1_ - it.it1_; 01255 } 01256 01257 // Dereference 01258 BOOST_UBLAS_INLINE 01259 const_reference operator * () const { 01260 size_type i = index1 (); 01261 size_type j = index2 (); 01262 BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ()); 01263 BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ()); 01264 if (triangular_type::other (i, j)) 01265 return *it1_; 01266 else 01267 return (*this) () (i, j); 01268 } 01269 BOOST_UBLAS_INLINE 01270 const_reference operator [] (difference_type n) const { 01271 return *(*this + n); 01272 } 01273 01274 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 01275 BOOST_UBLAS_INLINE 01276 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 01277 typename self_type:: 01278 #endif 01279 const_iterator2 begin () const { 01280 return (*this) ().find2 (1, index1 (), 0); 01281 } 01282 BOOST_UBLAS_INLINE 01283 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 01284 typename self_type:: 01285 #endif 01286 const_iterator2 end () const { 01287 return (*this) ().find2 (1, index1 (), (*this) ().size2 ()); 01288 } 01289 BOOST_UBLAS_INLINE 01290 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 01291 typename self_type:: 01292 #endif 01293 const_reverse_iterator2 rbegin () const { 01294 return const_reverse_iterator2 (end ()); 01295 } 01296 BOOST_UBLAS_INLINE 01297 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 01298 typename self_type:: 01299 #endif 01300 const_reverse_iterator2 rend () const { 01301 return const_reverse_iterator2 (begin ()); 01302 } 01303 #endif 01304 01305 // Indices 01306 BOOST_UBLAS_INLINE 01307 size_type index1 () const { 01308 return it1_.index1 (); 01309 } 01310 BOOST_UBLAS_INLINE 01311 size_type index2 () const { 01312 return it1_.index2 (); 01313 } 01314 01315 // Assignment 01316 BOOST_UBLAS_INLINE 01317 const_iterator1 &operator = (const const_iterator1 &it) { 01318 container_const_reference<self_type>::assign (&it ()); 01319 it1_ = it.it1_; 01320 return *this; 01321 } 01322 01323 // Comparison 01324 BOOST_UBLAS_INLINE 01325 bool operator == (const const_iterator1 &it) const { 01326 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 01327 return it1_ == it.it1_; 01328 } 01329 BOOST_UBLAS_INLINE 01330 bool operator < (const const_iterator1 &it) const { 01331 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 01332 return it1_ < it.it1_; 01333 } 01334 01335 private: 01336 const_subiterator1_type it1_; 01337 }; 01338 #endif 01339 01340 BOOST_UBLAS_INLINE 01341 const_iterator1 begin1 () const { 01342 return find1 (0, 0, 0); 01343 } 01344 BOOST_UBLAS_INLINE 01345 const_iterator1 end1 () const { 01346 return find1 (0, size1 (), 0); 01347 } 01348 01349 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR 01350 class iterator1: 01351 public container_reference<triangular_adaptor>, 01352 public random_access_iterator_base<typename iterator_restrict_traits< 01353 typename subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category, 01354 iterator1, value_type> { 01355 public: 01356 typedef typename subiterator1_type::value_type value_type; 01357 typedef typename subiterator1_type::difference_type difference_type; 01358 typedef typename subiterator1_type::reference reference; 01359 typedef typename subiterator1_type::pointer pointer; 01360 01361 typedef iterator2 dual_iterator_type; 01362 typedef reverse_iterator2 dual_reverse_iterator_type; 01363 01364 // Construction and destruction 01365 BOOST_UBLAS_INLINE 01366 iterator1 (): 01367 container_reference<self_type> (), it1_ () {} 01368 BOOST_UBLAS_INLINE 01369 iterator1 (self_type &m, const subiterator1_type &it1): 01370 container_reference<self_type> (m), it1_ (it1) {} 01371 01372 // Arithmetic 01373 BOOST_UBLAS_INLINE 01374 iterator1 &operator ++ () { 01375 ++ it1_; 01376 return *this; 01377 } 01378 BOOST_UBLAS_INLINE 01379 iterator1 &operator -- () { 01380 -- it1_; 01381 return *this; 01382 } 01383 BOOST_UBLAS_INLINE 01384 iterator1 &operator += (difference_type n) { 01385 it1_ += n; 01386 return *this; 01387 } 01388 BOOST_UBLAS_INLINE 01389 iterator1 &operator -= (difference_type n) { 01390 it1_ -= n; 01391 return *this; 01392 } 01393 BOOST_UBLAS_INLINE 01394 difference_type operator - (const iterator1 &it) const { 01395 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 01396 return it1_ - it.it1_; 01397 } 01398 01399 // Dereference 01400 BOOST_UBLAS_INLINE 01401 reference operator * () const { 01402 size_type i = index1 (); 01403 size_type j = index2 (); 01404 BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ()); 01405 BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ()); 01406 if (triangular_type::other (i, j)) 01407 return *it1_; 01408 else 01409 return (*this) () (i, j); 01410 } 01411 BOOST_UBLAS_INLINE 01412 reference operator [] (difference_type n) const { 01413 return *(*this + n); 01414 } 01415 01416 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 01417 BOOST_UBLAS_INLINE 01418 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 01419 typename self_type:: 01420 #endif 01421 iterator2 begin () const { 01422 return (*this) ().find2 (1, index1 (), 0); 01423 } 01424 BOOST_UBLAS_INLINE 01425 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 01426 typename self_type:: 01427 #endif 01428 iterator2 end () const { 01429 return (*this) ().find2 (1, index1 (), (*this) ().size2 ()); 01430 } 01431 BOOST_UBLAS_INLINE 01432 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 01433 typename self_type:: 01434 #endif 01435 reverse_iterator2 rbegin () const { 01436 return reverse_iterator2 (end ()); 01437 } 01438 BOOST_UBLAS_INLINE 01439 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 01440 typename self_type:: 01441 #endif 01442 reverse_iterator2 rend () const { 01443 return reverse_iterator2 (begin ()); 01444 } 01445 #endif 01446 01447 // Indices 01448 BOOST_UBLAS_INLINE 01449 size_type index1 () const { 01450 return it1_.index1 (); 01451 } 01452 BOOST_UBLAS_INLINE 01453 size_type index2 () const { 01454 return it1_.index2 (); 01455 } 01456 01457 // Assignment 01458 BOOST_UBLAS_INLINE 01459 iterator1 &operator = (const iterator1 &it) { 01460 container_reference<self_type>::assign (&it ()); 01461 it1_ = it.it1_; 01462 return *this; 01463 } 01464 01465 // Comparison 01466 BOOST_UBLAS_INLINE 01467 bool operator == (const iterator1 &it) const { 01468 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 01469 return it1_ == it.it1_; 01470 } 01471 BOOST_UBLAS_INLINE 01472 bool operator < (const iterator1 &it) const { 01473 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 01474 return it1_ < it.it1_; 01475 } 01476 01477 private: 01478 subiterator1_type it1_; 01479 01480 friend class const_iterator1; 01481 }; 01482 #endif 01483 01484 BOOST_UBLAS_INLINE 01485 iterator1 begin1 () { 01486 return find1 (0, 0, 0); 01487 } 01488 BOOST_UBLAS_INLINE 01489 iterator1 end1 () { 01490 return find1 (0, size1 (), 0); 01491 } 01492 01493 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR 01494 class const_iterator2: 01495 public container_const_reference<triangular_adaptor>, 01496 public random_access_iterator_base<typename iterator_restrict_traits< 01497 typename const_subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category, 01498 const_iterator2, value_type> { 01499 public: 01500 typedef typename const_subiterator2_type::value_type value_type; 01501 typedef typename const_subiterator2_type::difference_type difference_type; 01502 typedef typename const_subiterator2_type::reference reference; 01503 typedef typename const_subiterator2_type::pointer pointer; 01504 01505 typedef const_iterator1 dual_iterator_type; 01506 typedef const_reverse_iterator1 dual_reverse_iterator_type; 01507 01508 // Construction and destruction 01509 BOOST_UBLAS_INLINE 01510 const_iterator2 (): 01511 container_const_reference<self_type> (), it2_ () {} 01512 BOOST_UBLAS_INLINE 01513 const_iterator2 (const self_type &m, const const_subiterator2_type &it2): 01514 container_const_reference<self_type> (m), it2_ (it2) {} 01515 BOOST_UBLAS_INLINE 01516 const_iterator2 (const iterator2 &it): 01517 container_const_reference<self_type> (it ()), it2_ (it.it2_) {} 01518 01519 // Arithmetic 01520 BOOST_UBLAS_INLINE 01521 const_iterator2 &operator ++ () { 01522 ++ it2_; 01523 return *this; 01524 } 01525 BOOST_UBLAS_INLINE 01526 const_iterator2 &operator -- () { 01527 -- it2_; 01528 return *this; 01529 } 01530 BOOST_UBLAS_INLINE 01531 const_iterator2 &operator += (difference_type n) { 01532 it2_ += n; 01533 return *this; 01534 } 01535 BOOST_UBLAS_INLINE 01536 const_iterator2 &operator -= (difference_type n) { 01537 it2_ -= n; 01538 return *this; 01539 } 01540 BOOST_UBLAS_INLINE 01541 difference_type operator - (const const_iterator2 &it) const { 01542 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 01543 return it2_ - it.it2_; 01544 } 01545 01546 // Dereference 01547 BOOST_UBLAS_INLINE 01548 const_reference operator * () const { 01549 size_type i = index1 (); 01550 size_type j = index2 (); 01551 BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ()); 01552 BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ()); 01553 if (triangular_type::other (i, j)) 01554 return *it2_; 01555 else 01556 return (*this) () (i, j); 01557 } 01558 BOOST_UBLAS_INLINE 01559 const_reference operator [] (difference_type n) const { 01560 return *(*this + n); 01561 } 01562 01563 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 01564 BOOST_UBLAS_INLINE 01565 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 01566 typename self_type:: 01567 #endif 01568 const_iterator1 begin () const { 01569 return (*this) ().find1 (1, 0, index2 ()); 01570 } 01571 BOOST_UBLAS_INLINE 01572 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 01573 typename self_type:: 01574 #endif 01575 const_iterator1 end () const { 01576 return (*this) ().find1 (1, (*this) ().size1 (), index2 ()); 01577 } 01578 BOOST_UBLAS_INLINE 01579 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 01580 typename self_type:: 01581 #endif 01582 const_reverse_iterator1 rbegin () const { 01583 return const_reverse_iterator1 (end ()); 01584 } 01585 BOOST_UBLAS_INLINE 01586 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 01587 typename self_type:: 01588 #endif 01589 const_reverse_iterator1 rend () const { 01590 return const_reverse_iterator1 (begin ()); 01591 } 01592 #endif 01593 01594 // Indices 01595 BOOST_UBLAS_INLINE 01596 size_type index1 () const { 01597 return it2_.index1 (); 01598 } 01599 BOOST_UBLAS_INLINE 01600 size_type index2 () const { 01601 return it2_.index2 (); 01602 } 01603 01604 // Assignment 01605 BOOST_UBLAS_INLINE 01606 const_iterator2 &operator = (const const_iterator2 &it) { 01607 container_const_reference<self_type>::assign (&it ()); 01608 it2_ = it.it2_; 01609 return *this; 01610 } 01611 01612 // Comparison 01613 BOOST_UBLAS_INLINE 01614 bool operator == (const const_iterator2 &it) const { 01615 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 01616 return it2_ == it.it2_; 01617 } 01618 BOOST_UBLAS_INLINE 01619 bool operator < (const const_iterator2 &it) const { 01620 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 01621 return it2_ < it.it2_; 01622 } 01623 01624 private: 01625 const_subiterator2_type it2_; 01626 }; 01627 #endif 01628 01629 BOOST_UBLAS_INLINE 01630 const_iterator2 begin2 () const { 01631 return find2 (0, 0, 0); 01632 } 01633 BOOST_UBLAS_INLINE 01634 const_iterator2 end2 () const { 01635 return find2 (0, 0, size2 ()); 01636 } 01637 01638 #ifndef BOOST_UBLAS_USE_INDEXED_ITERATOR 01639 class iterator2: 01640 public container_reference<triangular_adaptor>, 01641 public random_access_iterator_base<typename iterator_restrict_traits< 01642 typename subiterator1_type::iterator_category, packed_random_access_iterator_tag>::iterator_category, 01643 iterator2, value_type> { 01644 public: 01645 typedef typename subiterator2_type::value_type value_type; 01646 typedef typename subiterator2_type::difference_type difference_type; 01647 typedef typename subiterator2_type::reference reference; 01648 typedef typename subiterator2_type::pointer pointer; 01649 01650 typedef iterator1 dual_iterator_type; 01651 typedef reverse_iterator1 dual_reverse_iterator_type; 01652 01653 // Construction and destruction 01654 BOOST_UBLAS_INLINE 01655 iterator2 (): 01656 container_reference<self_type> (), it2_ () {} 01657 BOOST_UBLAS_INLINE 01658 iterator2 (self_type &m, const subiterator2_type &it2): 01659 container_reference<self_type> (m), it2_ (it2) {} 01660 01661 // Arithmetic 01662 BOOST_UBLAS_INLINE 01663 iterator2 &operator ++ () { 01664 ++ it2_; 01665 return *this; 01666 } 01667 BOOST_UBLAS_INLINE 01668 iterator2 &operator -- () { 01669 -- it2_; 01670 return *this; 01671 } 01672 BOOST_UBLAS_INLINE 01673 iterator2 &operator += (difference_type n) { 01674 it2_ += n; 01675 return *this; 01676 } 01677 BOOST_UBLAS_INLINE 01678 iterator2 &operator -= (difference_type n) { 01679 it2_ -= n; 01680 return *this; 01681 } 01682 BOOST_UBLAS_INLINE 01683 difference_type operator - (const iterator2 &it) const { 01684 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 01685 return it2_ - it.it2_; 01686 } 01687 01688 // Dereference 01689 BOOST_UBLAS_INLINE 01690 reference operator * () const { 01691 size_type i = index1 (); 01692 size_type j = index2 (); 01693 BOOST_UBLAS_CHECK (i < (*this) ().size1 (), bad_index ()); 01694 BOOST_UBLAS_CHECK (j < (*this) ().size2 (), bad_index ()); 01695 if (triangular_type::other (i, j)) 01696 return *it2_; 01697 else 01698 return (*this) () (i, j); 01699 } 01700 BOOST_UBLAS_INLINE 01701 reference operator [] (difference_type n) const { 01702 return *(*this + n); 01703 } 01704 01705 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 01706 BOOST_UBLAS_INLINE 01707 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 01708 typename self_type:: 01709 #endif 01710 iterator1 begin () const { 01711 return (*this) ().find1 (1, 0, index2 ()); 01712 } 01713 BOOST_UBLAS_INLINE 01714 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 01715 typename self_type:: 01716 #endif 01717 iterator1 end () const { 01718 return (*this) ().find1 (1, (*this) ().size1 (), index2 ()); 01719 } 01720 BOOST_UBLAS_INLINE 01721 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 01722 typename self_type:: 01723 #endif 01724 reverse_iterator1 rbegin () const { 01725 return reverse_iterator1 (end ()); 01726 } 01727 BOOST_UBLAS_INLINE 01728 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 01729 typename self_type:: 01730 #endif 01731 reverse_iterator1 rend () const { 01732 return reverse_iterator1 (begin ()); 01733 } 01734 #endif 01735 01736 // Indices 01737 BOOST_UBLAS_INLINE 01738 size_type index1 () const { 01739 return it2_.index1 (); 01740 } 01741 BOOST_UBLAS_INLINE 01742 size_type index2 () const { 01743 return it2_.index2 (); 01744 } 01745 01746 // Assignment 01747 BOOST_UBLAS_INLINE 01748 iterator2 &operator = (const iterator2 &it) { 01749 container_reference<self_type>::assign (&it ()); 01750 it2_ = it.it2_; 01751 return *this; 01752 } 01753 01754 // Comparison 01755 BOOST_UBLAS_INLINE 01756 bool operator == (const iterator2 &it) const { 01757 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 01758 return it2_ == it.it2_; 01759 } 01760 BOOST_UBLAS_INLINE 01761 bool operator < (const iterator2 &it) const { 01762 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 01763 return it2_ < it.it2_; 01764 } 01765 01766 private: 01767 subiterator2_type it2_; 01768 01769 friend class const_iterator2; 01770 }; 01771 #endif 01772 01773 BOOST_UBLAS_INLINE 01774 iterator2 begin2 () { 01775 return find2 (0, 0, 0); 01776 } 01777 BOOST_UBLAS_INLINE 01778 iterator2 end2 () { 01779 return find2 (0, 0, size2 ()); 01780 } 01781 01782 // Reverse iterators 01783 01784 BOOST_UBLAS_INLINE 01785 const_reverse_iterator1 rbegin1 () const { 01786 return const_reverse_iterator1 (end1 ()); 01787 } 01788 BOOST_UBLAS_INLINE 01789 const_reverse_iterator1 rend1 () const { 01790 return const_reverse_iterator1 (begin1 ()); 01791 } 01792 01793 BOOST_UBLAS_INLINE 01794 reverse_iterator1 rbegin1 () { 01795 return reverse_iterator1 (end1 ()); 01796 } 01797 BOOST_UBLAS_INLINE 01798 reverse_iterator1 rend1 () { 01799 return reverse_iterator1 (begin1 ()); 01800 } 01801 01802 BOOST_UBLAS_INLINE 01803 const_reverse_iterator2 rbegin2 () const { 01804 return const_reverse_iterator2 (end2 ()); 01805 } 01806 BOOST_UBLAS_INLINE 01807 const_reverse_iterator2 rend2 () const { 01808 return const_reverse_iterator2 (begin2 ()); 01809 } 01810 01811 BOOST_UBLAS_INLINE 01812 reverse_iterator2 rbegin2 () { 01813 return reverse_iterator2 (end2 ()); 01814 } 01815 BOOST_UBLAS_INLINE 01816 reverse_iterator2 rend2 () { 01817 return reverse_iterator2 (begin2 ()); 01818 } 01819 01820 private: 01821 matrix_closure_type data_; 01822 static const value_type zero_; 01823 static const value_type one_; 01824 }; 01825 01826 template<class M, class TRI> 01827 const typename triangular_adaptor<M, TRI>::value_type triangular_adaptor<M, TRI>::zero_ = value_type/*zero*/(); 01828 template<class M, class TRI> 01829 const typename triangular_adaptor<M, TRI>::value_type triangular_adaptor<M, TRI>::one_ (1); 01830 01831 template <class M, class TRI> 01832 struct vector_temporary_traits< triangular_adaptor<M, TRI> > 01833 : vector_temporary_traits< typename boost::remove_const<M>::type > {} ; 01834 template <class M, class TRI> 01835 struct vector_temporary_traits< const triangular_adaptor<M, TRI> > 01836 : vector_temporary_traits< typename boost::remove_const<M>::type > {} ; 01837 01838 template <class M, class TRI> 01839 struct matrix_temporary_traits< triangular_adaptor<M, TRI> > 01840 : matrix_temporary_traits< typename boost::remove_const<M>::type > {}; 01841 template <class M, class TRI> 01842 struct matrix_temporary_traits< const triangular_adaptor<M, TRI> > 01843 : matrix_temporary_traits< typename boost::remove_const<M>::type > {}; 01844 01845 01846 template<class E1, class E2> 01847 struct matrix_vector_solve_traits { 01848 typedef typename promote_traits<typename E1::value_type, typename E2::value_type>::promote_type promote_type; 01849 typedef vector<promote_type> result_type; 01850 }; 01851 01852 // Operations: 01853 // n * (n - 1) / 2 + n = n * (n + 1) / 2 multiplications, 01854 // n * (n - 1) / 2 additions 01855 01856 // Dense (proxy) case 01857 template<class E1, class E2> 01858 BOOST_UBLAS_INLINE 01859 void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2, 01860 lower_tag, column_major_tag, dense_proxy_tag) { 01861 typedef typename E2::size_type size_type; 01862 typedef typename E2::difference_type difference_type; 01863 typedef typename E2::value_type value_type; 01864 01865 BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ()); 01866 BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ()); 01867 size_type size = e2 ().size (); 01868 for (size_type n = 0; n < size; ++ n) { 01869 #ifndef BOOST_UBLAS_SINGULAR_CHECK 01870 BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ()); 01871 #else 01872 if (e1 () (n, n) == value_type/*zero*/()) 01873 singular ().raise (); 01874 #endif 01875 value_type t = e2 () (n) /= e1 () (n, n); 01876 if (t != value_type/*zero*/()) { 01877 for (size_type m = n + 1; m < size; ++ m) 01878 e2 () (m) -= e1 () (m, n) * t; 01879 } 01880 } 01881 } 01882 // Packed (proxy) case 01883 template<class E1, class E2> 01884 BOOST_UBLAS_INLINE 01885 void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2, 01886 lower_tag, column_major_tag, packed_proxy_tag) { 01887 typedef typename E2::size_type size_type; 01888 typedef typename E2::difference_type difference_type; 01889 typedef typename E2::value_type value_type; 01890 01891 BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ()); 01892 BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ()); 01893 size_type size = e2 ().size (); 01894 for (size_type n = 0; n < size; ++ n) { 01895 #ifndef BOOST_UBLAS_SINGULAR_CHECK 01896 BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ()); 01897 #else 01898 if (e1 () (n, n) == value_type/*zero*/()) 01899 singular ().raise (); 01900 #endif 01901 value_type t = e2 () (n) /= e1 () (n, n); 01902 if (t != value_type/*zero*/()) { 01903 typename E1::const_iterator1 it1e1 (e1 ().find1 (1, n + 1, n)); 01904 typename E1::const_iterator1 it1e1_end (e1 ().find1 (1, e1 ().size1 (), n)); 01905 difference_type m (it1e1_end - it1e1); 01906 while (-- m >= 0) 01907 e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1; 01908 } 01909 } 01910 } 01911 // Sparse (proxy) case 01912 template<class E1, class E2> 01913 BOOST_UBLAS_INLINE 01914 void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2, 01915 lower_tag, column_major_tag, unknown_storage_tag) { 01916 typedef typename E2::size_type size_type; 01917 typedef typename E2::difference_type difference_type; 01918 typedef typename E2::value_type value_type; 01919 01920 BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ()); 01921 BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ()); 01922 size_type size = e2 ().size (); 01923 for (size_type n = 0; n < size; ++ n) { 01924 #ifndef BOOST_UBLAS_SINGULAR_CHECK 01925 BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ()); 01926 #else 01927 if (e1 () (n, n) == value_type/*zero*/()) 01928 singular ().raise (); 01929 #endif 01930 value_type t = e2 () (n) /= e1 () (n, n); 01931 if (t != value_type/*zero*/()) { 01932 typename E1::const_iterator1 it1e1 (e1 ().find1 (1, n + 1, n)); 01933 typename E1::const_iterator1 it1e1_end (e1 ().find1 (1, e1 ().size1 (), n)); 01934 while (it1e1 != it1e1_end) 01935 e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1; 01936 } 01937 } 01938 } 01939 // Redirectors :-) 01940 template<class E1, class E2> 01941 BOOST_UBLAS_INLINE 01942 void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2, 01943 lower_tag, column_major_tag) { 01944 typedef typename E1::storage_category storage_category; 01945 inplace_solve (e1, e2, 01946 lower_tag (), column_major_tag (), storage_category ()); 01947 } 01948 template<class E1, class E2> 01949 BOOST_UBLAS_INLINE 01950 void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2, 01951 lower_tag, row_major_tag) { 01952 typedef typename E1::storage_category storage_category; 01953 inplace_solve (e2, trans (e1), 01954 upper_tag (), row_major_tag (), storage_category ()); 01955 } 01956 // Dispatcher 01957 template<class E1, class E2> 01958 BOOST_UBLAS_INLINE 01959 void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2, 01960 lower_tag) { 01961 typedef typename E1::orientation_category orientation_category; 01962 inplace_solve (e1, e2, 01963 lower_tag (), orientation_category ()); 01964 } 01965 template<class E1, class E2> 01966 BOOST_UBLAS_INLINE 01967 void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2, 01968 unit_lower_tag) { 01969 typedef typename E1::orientation_category orientation_category; 01970 inplace_solve (triangular_adaptor<const E1, unit_lower> (e1 ()), e2, 01971 unit_lower_tag (), orientation_category ()); 01972 } 01973 01974 // Dense (proxy) case 01975 template<class E1, class E2> 01976 BOOST_UBLAS_INLINE 01977 void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2, 01978 upper_tag, column_major_tag, dense_proxy_tag) { 01979 typedef typename E2::size_type size_type; 01980 typedef typename E2::difference_type difference_type; 01981 typedef typename E2::value_type value_type; 01982 01983 BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ()); 01984 BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ()); 01985 size_type size = e2 ().size (); 01986 for (difference_type n = size - 1; n >= 0; -- n) { 01987 #ifndef BOOST_UBLAS_SINGULAR_CHECK 01988 BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ()); 01989 #else 01990 if (e1 () (n, n) == value_type/*zero*/()) 01991 singular ().raise (); 01992 #endif 01993 value_type t = e2 () (n) /= e1 () (n, n); 01994 if (t != value_type/*zero*/()) { 01995 for (difference_type m = n - 1; m >= 0; -- m) 01996 e2 () (m) -= e1 () (m, n) * t; 01997 } 01998 } 01999 } 02000 // Packed (proxy) case 02001 template<class E1, class E2> 02002 BOOST_UBLAS_INLINE 02003 void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2, 02004 upper_tag, column_major_tag, packed_proxy_tag) { 02005 typedef typename E2::size_type size_type; 02006 typedef typename E2::difference_type difference_type; 02007 typedef typename E2::value_type value_type; 02008 02009 BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ()); 02010 BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ()); 02011 size_type size = e2 ().size (); 02012 for (difference_type n = size - 1; n >= 0; -- n) { 02013 #ifndef BOOST_UBLAS_SINGULAR_CHECK 02014 BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ()); 02015 #else 02016 if (e1 () (n, n) == value_type/*zero*/()) 02017 singular ().raise (); 02018 #endif 02019 value_type t = e2 () (n) /= e1 () (n, n); 02020 if (t != value_type/*zero*/()) { 02021 typename E1::const_reverse_iterator1 it1e1 (e1 ().find1 (1, n, n)); 02022 typename E1::const_reverse_iterator1 it1e1_rend (e1 ().find1 (1, 0, n)); 02023 difference_type m (it1e1_rend - it1e1); 02024 while (-- m >= 0) 02025 e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1; 02026 } 02027 } 02028 } 02029 // Sparse (proxy) case 02030 template<class E1, class E2> 02031 BOOST_UBLAS_INLINE 02032 void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2, 02033 upper_tag, column_major_tag, unknown_storage_tag) { 02034 typedef typename E2::size_type size_type; 02035 typedef typename E2::difference_type difference_type; 02036 typedef typename E2::value_type value_type; 02037 02038 BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ()); 02039 BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size (), bad_size ()); 02040 size_type size = e2 ().size (); 02041 for (difference_type n = size - 1; n >= 0; -- n) { 02042 #ifndef BOOST_UBLAS_SINGULAR_CHECK 02043 BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ()); 02044 #else 02045 if (e1 () (n, n) == value_type/*zero*/()) 02046 singular ().raise (); 02047 #endif 02048 value_type t = e2 () (n) /= e1 () (n, n); 02049 if (t != value_type/*zero*/()) { 02050 typename E1::const_reverse_iterator1 it1e1 (e1 ().find1 (1, n, n)); 02051 typename E1::const_reverse_iterator1 it1e1_rend (e1 ().find1 (1, 0, n)); 02052 while (it1e1 != it1e1_rend) 02053 e2 () (it1e1.index1 ()) -= *it1e1 * t, ++ it1e1; 02054 } 02055 } 02056 } 02057 // Redirectors :-) 02058 template<class E1, class E2> 02059 BOOST_UBLAS_INLINE 02060 void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2, 02061 upper_tag, column_major_tag) { 02062 typedef typename E1::storage_category storage_category; 02063 inplace_solve (e1, e2, 02064 upper_tag (), column_major_tag (), storage_category ()); 02065 } 02066 template<class E1, class E2> 02067 BOOST_UBLAS_INLINE 02068 void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2, 02069 upper_tag, row_major_tag) { 02070 typedef typename E1::storage_category storage_category; 02071 inplace_solve (e2, trans (e1), 02072 lower_tag (), row_major_tag (), storage_category ()); 02073 } 02074 // Dispatcher 02075 template<class E1, class E2> 02076 BOOST_UBLAS_INLINE 02077 void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2, 02078 upper_tag) { 02079 typedef typename E1::orientation_category orientation_category; 02080 inplace_solve (e1, e2, 02081 upper_tag (), orientation_category ()); 02082 } 02083 template<class E1, class E2> 02084 BOOST_UBLAS_INLINE 02085 void inplace_solve (const matrix_expression<E1> &e1, vector_expression<E2> &e2, 02086 unit_upper_tag) { 02087 typedef typename E1::orientation_category orientation_category; 02088 inplace_solve (triangular_adaptor<const E1, unit_upper> (e1 ()), e2, 02089 unit_upper_tag (), orientation_category ()); 02090 } 02091 02092 template<class E1, class E2, class C> 02093 BOOST_UBLAS_INLINE 02094 typename matrix_vector_solve_traits<E1, E2>::result_type 02095 solve (const matrix_expression<E1> &e1, 02096 const vector_expression<E2> &e2, 02097 C) { 02098 typename matrix_vector_solve_traits<E1, E2>::result_type r (e2); 02099 inplace_solve (e1, r, C ()); 02100 return r; 02101 } 02102 02103 // Dense (proxy) case 02104 template<class E1, class E2> 02105 BOOST_UBLAS_INLINE 02106 void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2, 02107 lower_tag, row_major_tag, dense_proxy_tag) { 02108 typedef typename E1::size_type size_type; 02109 typedef typename E1::difference_type difference_type; 02110 typedef typename E1::value_type value_type; 02111 02112 BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ()); 02113 BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ()); 02114 size_type size = e1 ().size (); 02115 for (difference_type n = size - 1; n >= 0; -- n) { 02116 #ifndef BOOST_UBLAS_SINGULAR_CHECK 02117 BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ()); 02118 #else 02119 if (e2 () (n, n) == value_type/*zero*/()) 02120 singular ().raise (); 02121 #endif 02122 value_type t = e1 () (n) /= e2 () (n, n); 02123 if (t != value_type/*zero*/()) { 02124 for (difference_type m = n - 1; m >= 0; -- m) 02125 e1 () (m) -= t * e2 () (n, m); 02126 } 02127 } 02128 } 02129 // Packed (proxy) case 02130 template<class E1, class E2> 02131 BOOST_UBLAS_INLINE 02132 void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2, 02133 lower_tag, row_major_tag, packed_proxy_tag) { 02134 typedef typename E1::size_type size_type; 02135 typedef typename E1::difference_type difference_type; 02136 typedef typename E1::value_type value_type; 02137 02138 BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ()); 02139 BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ()); 02140 size_type size = e1 ().size (); 02141 for (difference_type n = size - 1; n >= 0; -- n) { 02142 #ifndef BOOST_UBLAS_SINGULAR_CHECK 02143 BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ()); 02144 #else 02145 if (e2 () (n, n) == value_type/*zero*/()) 02146 singular ().raise (); 02147 #endif 02148 value_type t = e1 () (n) /= e2 () (n, n); 02149 if (t != value_type/*zero*/()) { 02150 typename E2::const_reverse_iterator2 it2e2 (e2 ().find2 (1, n, n)); 02151 typename E2::const_reverse_iterator2 it2e2_rend (e2 ().find2 (1, n, 0)); 02152 difference_type m (it2e2_rend - it2e2); 02153 while (-- m >= 0) 02154 e1 () (it2e2.index2 ()) -= *it2e2 * t, ++ it2e2; 02155 } 02156 } 02157 } 02158 // Sparse (proxy) case 02159 template<class E1, class E2> 02160 BOOST_UBLAS_INLINE 02161 void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2, 02162 lower_tag, row_major_tag, unknown_storage_tag) { 02163 typedef typename E1::size_type size_type; 02164 typedef typename E1::difference_type difference_type; 02165 typedef typename E1::value_type value_type; 02166 02167 BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ()); 02168 BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ()); 02169 size_type size = e1 ().size (); 02170 for (difference_type n = size - 1; n >= 0; -- n) { 02171 #ifndef BOOST_UBLAS_SINGULAR_CHECK 02172 BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ()); 02173 #else 02174 if (e2 () (n, n) == value_type/*zero*/()) 02175 singular ().raise (); 02176 #endif 02177 value_type t = e1 () (n) /= e2 () (n, n); 02178 if (t != value_type/*zero*/()) { 02179 typename E2::const_reverse_iterator2 it2e2 (e2 ().find2 (1, n, n)); 02180 typename E2::const_reverse_iterator2 it2e2_rend (e2 ().find2 (1, n, 0)); 02181 while (it2e2 != it2e2_rend) 02182 e1 () (it2e2.index2 ()) -= *it2e2 * t, ++ it2e2; 02183 } 02184 } 02185 } 02186 // Redirectors :-) 02187 template<class E1, class E2> 02188 BOOST_UBLAS_INLINE 02189 void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2, 02190 lower_tag, row_major_tag) { 02191 typedef typename E1::storage_category storage_category; 02192 inplace_solve (e1, e2, 02193 lower_tag (), row_major_tag (), storage_category ()); 02194 } 02195 template<class E1, class E2> 02196 BOOST_UBLAS_INLINE 02197 void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2, 02198 lower_tag, column_major_tag) { 02199 typedef typename E1::storage_category storage_category; 02200 inplace_solve (trans (e2), e1, 02201 upper_tag (), row_major_tag (), storage_category ()); 02202 } 02203 // Dispatcher 02204 template<class E1, class E2> 02205 BOOST_UBLAS_INLINE 02206 void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2, 02207 lower_tag) { 02208 typedef typename E2::orientation_category orientation_category; 02209 inplace_solve (e1, e2, 02210 lower_tag (), orientation_category ()); 02211 } 02212 template<class E1, class E2> 02213 BOOST_UBLAS_INLINE 02214 void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2, 02215 unit_lower_tag) { 02216 typedef typename E2::orientation_category orientation_category; 02217 inplace_solve (e1, triangular_adaptor<const E2, unit_lower> (e2 ()), 02218 unit_lower_tag (), orientation_category ()); 02219 } 02220 02221 // Dense (proxy) case 02222 template<class E1, class E2> 02223 BOOST_UBLAS_INLINE 02224 void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2, 02225 upper_tag, row_major_tag, dense_proxy_tag) { 02226 typedef typename E1::size_type size_type; 02227 typedef typename E1::difference_type difference_type; 02228 typedef typename E1::value_type value_type; 02229 02230 BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ()); 02231 BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ()); 02232 size_type size = e1 ().size (); 02233 for (size_type n = 0; n < size; ++ n) { 02234 #ifndef BOOST_UBLAS_SINGULAR_CHECK 02235 BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ()); 02236 #else 02237 if (e2 () (n, n) == value_type/*zero*/()) 02238 singular ().raise (); 02239 #endif 02240 value_type t = e1 () (n) /= e2 () (n, n); 02241 if (t != value_type/*zero*/()) { 02242 for (size_type m = n + 1; m < size; ++ m) 02243 e1 () (m) -= t * e2 () (n, m); 02244 } 02245 } 02246 } 02247 // Packed (proxy) case 02248 template<class E1, class E2> 02249 BOOST_UBLAS_INLINE 02250 void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2, 02251 upper_tag, row_major_tag, packed_proxy_tag) { 02252 typedef typename E1::size_type size_type; 02253 typedef typename E1::difference_type difference_type; 02254 typedef typename E1::value_type value_type; 02255 02256 BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ()); 02257 BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ()); 02258 size_type size = e1 ().size (); 02259 for (size_type n = 0; n < size; ++ n) { 02260 #ifndef BOOST_UBLAS_SINGULAR_CHECK 02261 BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ()); 02262 #else 02263 if (e2 () (n, n) == value_type/*zero*/()) 02264 singular ().raise (); 02265 #endif 02266 value_type t = e1 () (n) /= e2 () (n, n); 02267 if (t != value_type/*zero*/()) { 02268 typename E2::const_iterator2 it2e2 (e2 ().find2 (1, n, n + 1)); 02269 typename E2::const_iterator2 it2e2_end (e2 ().find2 (1, n, e2 ().size2 ())); 02270 difference_type m (it2e2_end - it2e2); 02271 while (-- m >= 0) 02272 e1 () (it2e2.index2 ()) -= *it2e2 * t, ++ it2e2; 02273 } 02274 } 02275 } 02276 // Sparse (proxy) case 02277 template<class E1, class E2> 02278 BOOST_UBLAS_INLINE 02279 void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2, 02280 upper_tag, row_major_tag, unknown_storage_tag) { 02281 typedef typename E1::size_type size_type; 02282 typedef typename E1::difference_type difference_type; 02283 typedef typename E1::value_type value_type; 02284 02285 BOOST_UBLAS_CHECK (e1 ().size () == e2 ().size1 (), bad_size ()); 02286 BOOST_UBLAS_CHECK (e2 ().size1 () == e2 ().size2 (), bad_size ()); 02287 size_type size = e1 ().size (); 02288 for (size_type n = 0; n < size; ++ n) { 02289 #ifndef BOOST_UBLAS_SINGULAR_CHECK 02290 BOOST_UBLAS_CHECK (e2 () (n, n) != value_type/*zero*/(), singular ()); 02291 #else 02292 if (e2 () (n, n) == value_type/*zero*/()) 02293 singular ().raise (); 02294 #endif 02295 value_type t = e1 () (n) /= e2 () (n, n); 02296 if (t != value_type/*zero*/()) { 02297 typename E2::const_iterator2 it2e2 (e2 ().find2 (1, n, n + 1)); 02298 typename E2::const_iterator2 it2e2_end (e2 ().find2 (1, n, e2 ().size2 ())); 02299 while (it2e2 != it2e2_end) 02300 e1 () (it2e2.index2 ()) -= *it2e2 * t, ++ it2e2; 02301 } 02302 } 02303 } 02304 // Redirectors :-) 02305 template<class E1, class E2> 02306 BOOST_UBLAS_INLINE 02307 void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2, 02308 upper_tag, row_major_tag) { 02309 typedef typename E1::storage_category storage_category; 02310 inplace_solve (e1, e2, 02311 upper_tag (), row_major_tag (), storage_category ()); 02312 } 02313 template<class E1, class E2> 02314 BOOST_UBLAS_INLINE 02315 void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2, 02316 upper_tag, column_major_tag) { 02317 typedef typename E1::storage_category storage_category; 02318 inplace_solve (trans (e2), e1, 02319 lower_tag (), row_major_tag (), storage_category ()); 02320 } 02321 // Dispatcher 02322 template<class E1, class E2> 02323 BOOST_UBLAS_INLINE 02324 void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2, 02325 upper_tag) { 02326 typedef typename E2::orientation_category orientation_category; 02327 inplace_solve (e1, e2, 02328 upper_tag (), orientation_category ()); 02329 } 02330 template<class E1, class E2> 02331 BOOST_UBLAS_INLINE 02332 void inplace_solve (vector_expression<E1> &e1, const matrix_expression<E2> &e2, 02333 unit_upper_tag) { 02334 typedef typename E2::orientation_category orientation_category; 02335 inplace_solve (e1, triangular_adaptor<const E2, unit_upper> (e2 ()), 02336 unit_upper_tag (), orientation_category ()); 02337 } 02338 02339 template<class E1, class E2, class C> 02340 BOOST_UBLAS_INLINE 02341 typename matrix_vector_solve_traits<E1, E2>::result_type 02342 solve (const vector_expression<E1> &e1, 02343 const matrix_expression<E2> &e2, 02344 C) { 02345 typename matrix_vector_solve_traits<E1, E2>::result_type r (e1); 02346 inplace_solve (r, e2, C ()); 02347 return r; 02348 } 02349 02350 template<class E1, class E2> 02351 struct matrix_matrix_solve_traits { 02352 typedef typename promote_traits<typename E1::value_type, typename E2::value_type>::promote_type promote_type; 02353 typedef matrix<promote_type> result_type; 02354 }; 02355 02356 // Operations: 02357 // k * n * (n - 1) / 2 + k * n = k * n * (n + 1) / 2 multiplications, 02358 // k * n * (n - 1) / 2 additions 02359 02360 // Dense (proxy) case 02361 template<class E1, class E2> 02362 BOOST_UBLAS_INLINE 02363 void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2, 02364 lower_tag, dense_proxy_tag) { 02365 typedef typename E2::size_type size_type; 02366 typedef typename E2::difference_type difference_type; 02367 typedef typename E2::value_type value_type; 02368 02369 BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ()); 02370 BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ()); 02371 size_type size1 = e2 ().size1 (); 02372 size_type size2 = e2 ().size2 (); 02373 for (size_type n = 0; n < size1; ++ n) { 02374 #ifndef BOOST_UBLAS_SINGULAR_CHECK 02375 BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ()); 02376 #else 02377 if (e1 () (n, n) == value_type/*zero*/()) 02378 singular ().raise (); 02379 #endif 02380 for (size_type l = 0; l < size2; ++ l) { 02381 value_type t = e2 () (n, l) /= e1 () (n, n); 02382 if (t != value_type/*zero*/()) { 02383 for (size_type m = n + 1; m < size1; ++ m) 02384 e2 () (m, l) -= e1 () (m, n) * t; 02385 } 02386 } 02387 } 02388 } 02389 // Packed (proxy) case 02390 template<class E1, class E2> 02391 BOOST_UBLAS_INLINE 02392 void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2, 02393 lower_tag, packed_proxy_tag) { 02394 typedef typename E2::size_type size_type; 02395 typedef typename E2::difference_type difference_type; 02396 typedef typename E2::value_type value_type; 02397 02398 BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ()); 02399 BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ()); 02400 size_type size1 = e2 ().size1 (); 02401 size_type size2 = e2 ().size2 (); 02402 for (size_type n = 0; n < size1; ++ n) { 02403 #ifndef BOOST_UBLAS_SINGULAR_CHECK 02404 BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ()); 02405 #else 02406 if (e1 () (n, n) == value_type/*zero*/()) 02407 singular ().raise (); 02408 #endif 02409 for (size_type l = 0; l < size2; ++ l) { 02410 value_type t = e2 () (n, l) /= e1 () (n, n); 02411 if (t != value_type/*zero*/()) { 02412 typename E1::const_iterator1 it1e1 (e1 ().find1 (1, n + 1, n)); 02413 typename E1::const_iterator1 it1e1_end (e1 ().find1 (1, e1 ().size1 (), n)); 02414 difference_type m (it1e1_end - it1e1); 02415 while (-- m >= 0) 02416 e2 () (it1e1.index1 (), l) -= *it1e1 * t, ++ it1e1; 02417 } 02418 } 02419 } 02420 } 02421 // Sparse (proxy) case 02422 template<class E1, class E2> 02423 BOOST_UBLAS_INLINE 02424 void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2, 02425 lower_tag, unknown_storage_tag) { 02426 typedef typename E2::size_type size_type; 02427 typedef typename E2::difference_type difference_type; 02428 typedef typename E2::value_type value_type; 02429 02430 BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ()); 02431 BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ()); 02432 size_type size1 = e2 ().size1 (); 02433 size_type size2 = e2 ().size2 (); 02434 for (size_type n = 0; n < size1; ++ n) { 02435 #ifndef BOOST_UBLAS_SINGULAR_CHECK 02436 BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ()); 02437 #else 02438 if (e1 () (n, n) == value_type/*zero*/()) 02439 singular ().raise (); 02440 #endif 02441 for (size_type l = 0; l < size2; ++ l) { 02442 value_type t = e2 () (n, l) /= e1 () (n, n); 02443 if (t != value_type/*zero*/()) { 02444 typename E1::const_iterator1 it1e1 (e1 ().find1 (1, n + 1, n)); 02445 typename E1::const_iterator1 it1e1_end (e1 ().find1 (1, e1 ().size1 (), n)); 02446 while (it1e1 != it1e1_end) 02447 e2 () (it1e1.index1 (), l) -= *it1e1 * t, ++ it1e1; 02448 } 02449 } 02450 } 02451 } 02452 // Dispatcher 02453 template<class E1, class E2> 02454 BOOST_UBLAS_INLINE 02455 void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2, 02456 lower_tag) { 02457 typedef typename E1::storage_category dispatch_category; 02458 inplace_solve (e1, e2, 02459 lower_tag (), dispatch_category ()); 02460 } 02461 template<class E1, class E2> 02462 BOOST_UBLAS_INLINE 02463 void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2, 02464 unit_lower_tag) { 02465 typedef typename E1::storage_category dispatch_category; 02466 inplace_solve (triangular_adaptor<const E1, unit_lower> (e1 ()), e2, 02467 unit_lower_tag (), dispatch_category ()); 02468 } 02469 02470 // Dense (proxy) case 02471 template<class E1, class E2> 02472 BOOST_UBLAS_INLINE 02473 void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2, 02474 upper_tag, dense_proxy_tag) { 02475 typedef typename E2::size_type size_type; 02476 typedef typename E2::difference_type difference_type; 02477 typedef typename E2::value_type value_type; 02478 02479 BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ()); 02480 BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ()); 02481 size_type size1 = e2 ().size1 (); 02482 size_type size2 = e2 ().size2 (); 02483 for (difference_type n = size1 - 1; n >= 0; -- n) { 02484 #ifndef BOOST_UBLAS_SINGULAR_CHECK 02485 BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ()); 02486 #else 02487 if (e1 () (n, n) == value_type/*zero*/()) 02488 singular ().raise (); 02489 #endif 02490 for (difference_type l = size2 - 1; l >= 0; -- l) { 02491 value_type t = e2 () (n, l) /= e1 () (n, n); 02492 if (t != value_type/*zero*/()) { 02493 for (difference_type m = n - 1; m >= 0; -- m) 02494 e2 () (m, l) -= e1 () (m, n) * t; 02495 } 02496 } 02497 } 02498 } 02499 // Packed (proxy) case 02500 template<class E1, class E2> 02501 BOOST_UBLAS_INLINE 02502 void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2, 02503 upper_tag, packed_proxy_tag) { 02504 typedef typename E2::size_type size_type; 02505 typedef typename E2::difference_type difference_type; 02506 typedef typename E2::value_type value_type; 02507 02508 BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ()); 02509 BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ()); 02510 size_type size1 = e2 ().size1 (); 02511 size_type size2 = e2 ().size2 (); 02512 for (difference_type n = size1 - 1; n >= 0; -- n) { 02513 #ifndef BOOST_UBLAS_SINGULAR_CHECK 02514 BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ()); 02515 #else 02516 if (e1 () (n, n) == value_type/*zero*/()) 02517 singular ().raise (); 02518 #endif 02519 for (difference_type l = size2 - 1; l >= 0; -- l) { 02520 value_type t = e2 () (n, l) /= e1 () (n, n); 02521 if (t != value_type/*zero*/()) { 02522 typename E1::const_reverse_iterator1 it1e1 (e1 ().find1 (1, n, n)); 02523 typename E1::const_reverse_iterator1 it1e1_rend (e1 ().find1 (1, 0, n)); 02524 difference_type m (it1e1_rend - it1e1); 02525 while (-- m >= 0) 02526 e2 () (it1e1.index1 (), l) -= *it1e1 * t, ++ it1e1; 02527 } 02528 } 02529 } 02530 } 02531 // Sparse (proxy) case 02532 template<class E1, class E2> 02533 BOOST_UBLAS_INLINE 02534 void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2, 02535 upper_tag, unknown_storage_tag) { 02536 typedef typename E2::size_type size_type; 02537 typedef typename E2::difference_type difference_type; 02538 typedef typename E2::value_type value_type; 02539 02540 BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ()); 02541 BOOST_UBLAS_CHECK (e1 ().size2 () == e2 ().size1 (), bad_size ()); 02542 size_type size1 = e2 ().size1 (); 02543 size_type size2 = e2 ().size2 (); 02544 for (difference_type n = size1 - 1; n >= 0; -- n) { 02545 #ifndef BOOST_UBLAS_SINGULAR_CHECK 02546 BOOST_UBLAS_CHECK (e1 () (n, n) != value_type/*zero*/(), singular ()); 02547 #else 02548 if (e1 () (n, n) == value_type/*zero*/()) 02549 singular ().raise (); 02550 #endif 02551 for (difference_type l = size2 - 1; l >= 0; -- l) { 02552 value_type t = e2 () (n, l) /= e1 () (n, n); 02553 if (t != value_type/*zero*/()) { 02554 typename E1::const_reverse_iterator1 it1e1 (e1 ().find1 (1, n, n)); 02555 typename E1::const_reverse_iterator1 it1e1_rend (e1 ().find1 (1, 0, n)); 02556 while (it1e1 != it1e1_rend) 02557 e2 () (it1e1.index1 (), l) -= *it1e1 * t, ++ it1e1; 02558 } 02559 } 02560 } 02561 } 02562 // Dispatcher 02563 template<class E1, class E2> 02564 BOOST_UBLAS_INLINE 02565 void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2, 02566 upper_tag) { 02567 typedef typename E1::storage_category dispatch_category; 02568 inplace_solve (e1, e2, 02569 upper_tag (), dispatch_category ()); 02570 } 02571 template<class E1, class E2> 02572 BOOST_UBLAS_INLINE 02573 void inplace_solve (const matrix_expression<E1> &e1, matrix_expression<E2> &e2, 02574 unit_upper_tag) { 02575 typedef typename E1::storage_category dispatch_category; 02576 inplace_solve (triangular_adaptor<const E1, unit_upper> (e1 ()), e2, 02577 unit_upper_tag (), dispatch_category ()); 02578 } 02579 02580 template<class E1, class E2, class C> 02581 BOOST_UBLAS_INLINE 02582 typename matrix_matrix_solve_traits<E1, E2>::result_type 02583 solve (const matrix_expression<E1> &e1, 02584 const matrix_expression<E2> &e2, 02585 C) { 02586 typename matrix_matrix_solve_traits<E1, E2>::result_type r (e2); 02587 inplace_solve (e1, r, C ()); 02588 return r; 02589 } 02590 02591 }}} 02592 02593 #endif