1 // 2 // Copyright (c) 2000-2007 3 // Joerg Walter, Mathias Koch, Gunter Winkler 4 // 5 // Distributed under the Boost Software License, Version 1.0. (See 6 // accompanying file LICENSE_1_0.txt or copy at 7 // http://www.boost.org/LICENSE_1_0.txt) 8 // 9 // The authors gratefully acknowledge the support of 10 // GeNeSys mbH & Co. KG in producing this work. 11 // 12 13 #ifndef _BOOST_UBLAS_MATRIX_SPARSE_ 14 #define _BOOST_UBLAS_MATRIX_SPARSE_ 15 16 #include <boost/numeric/ublas/vector_sparse.hpp> 17 #include <boost/numeric/ublas/matrix_expression.hpp> 18 #include <boost/numeric/ublas/detail/matrix_assign.hpp> 19 #if BOOST_UBLAS_TYPE_CHECK 20 #include <boost/numeric/ublas/matrix.hpp> 21 #endif 22 23 // Iterators based on ideas of Jeremy Siek 24 25 namespace boost { namespace numeric { namespace ublas { 26 27 #ifdef BOOST_UBLAS_STRICT_MATRIX_SPARSE 28 29 template<class M> 30 class sparse_matrix_element: 31 public container_reference<M> { 32 public: 33 typedef M matrix_type; 34 typedef typename M::size_type size_type; 35 typedef typename M::value_type value_type; 36 typedef const value_type &const_reference; 37 typedef value_type *pointer; 38 typedef const value_type *const_pointer; 39 40 private: 41 // Proxied element operations get_d() const42 void get_d () const { 43 const_pointer p = (*this) ().find_element (i_, j_); 44 if (p) 45 d_ = *p; 46 else 47 d_ = value_type/*zero*/(); 48 } 49 set(const value_type & s) const50 void set (const value_type &s) const { 51 pointer p = (*this) ().find_element (i_, j_); 52 if (!p) 53 (*this) ().insert_element (i_, j_, s); 54 else 55 *p = s; 56 } 57 58 public: 59 // Construction and destruction 60 BOOST_UBLAS_INLINE sparse_matrix_element(matrix_type & m,size_type i,size_type j)61 sparse_matrix_element (matrix_type &m, size_type i, size_type j): 62 container_reference<matrix_type> (m), i_ (i), j_ (j) { 63 } 64 BOOST_UBLAS_INLINE sparse_matrix_element(const sparse_matrix_element & p)65 sparse_matrix_element (const sparse_matrix_element &p): 66 container_reference<matrix_type> (p), i_ (p.i_), j_ (p.j_) {} 67 BOOST_UBLAS_INLINE ~sparse_matrix_element()68 ~sparse_matrix_element () { 69 } 70 71 // Assignment 72 BOOST_UBLAS_INLINE operator =(const sparse_matrix_element & p)73 sparse_matrix_element &operator = (const sparse_matrix_element &p) { 74 // Overide the implict copy assignment 75 p.get_d (); 76 set (p.d_); 77 return *this; 78 } 79 template<class D> 80 BOOST_UBLAS_INLINE operator =(const D & d)81 sparse_matrix_element &operator = (const D &d) { 82 set (d); 83 return *this; 84 } 85 template<class D> 86 BOOST_UBLAS_INLINE operator +=(const D & d)87 sparse_matrix_element &operator += (const D &d) { 88 get_d (); 89 d_ += d; 90 set (d_); 91 return *this; 92 } 93 template<class D> 94 BOOST_UBLAS_INLINE operator -=(const D & d)95 sparse_matrix_element &operator -= (const D &d) { 96 get_d (); 97 d_ -= d; 98 set (d_); 99 return *this; 100 } 101 template<class D> 102 BOOST_UBLAS_INLINE operator *=(const D & d)103 sparse_matrix_element &operator *= (const D &d) { 104 get_d (); 105 d_ *= d; 106 set (d_); 107 return *this; 108 } 109 template<class D> 110 BOOST_UBLAS_INLINE operator /=(const D & d)111 sparse_matrix_element &operator /= (const D &d) { 112 get_d (); 113 d_ /= d; 114 set (d_); 115 return *this; 116 } 117 118 // Comparison 119 template<class D> 120 BOOST_UBLAS_INLINE operator ==(const D & d) const121 bool operator == (const D &d) const { 122 get_d (); 123 return d_ == d; 124 } 125 template<class D> 126 BOOST_UBLAS_INLINE operator !=(const D & d) const127 bool operator != (const D &d) const { 128 get_d (); 129 return d_ != d; 130 } 131 132 // Conversion - weak link in proxy as d_ is not a perfect alias for the element 133 BOOST_UBLAS_INLINE operator const_reference() const134 operator const_reference () const { 135 get_d (); 136 return d_; 137 } 138 139 // Conversion to reference - may be invalidated 140 BOOST_UBLAS_INLINE ref() const141 value_type& ref () const { 142 const pointer p = (*this) ().find_element (i_, j_); 143 if (!p) 144 return (*this) ().insert_element (i_, j_, value_type/*zero*/()); 145 else 146 return *p; 147 } 148 149 private: 150 size_type i_; 151 size_type j_; 152 mutable value_type d_; 153 }; 154 155 /* 156 * Generalise explicit reference access 157 */ 158 namespace detail { 159 template <class V> 160 struct element_reference<sparse_matrix_element<V> > { 161 typedef typename V::value_type& reference; get_referenceboost::numeric::ublas::detail::element_reference162 static reference get_reference (const sparse_matrix_element<V>& sve) 163 { 164 return sve.ref (); 165 } 166 }; 167 } 168 169 170 template<class M> 171 struct type_traits<sparse_matrix_element<M> > { 172 typedef typename M::value_type element_type; 173 typedef type_traits<sparse_matrix_element<M> > self_type; 174 typedef typename type_traits<element_type>::value_type value_type; 175 typedef typename type_traits<element_type>::const_reference const_reference; 176 typedef sparse_matrix_element<M> reference; 177 typedef typename type_traits<element_type>::real_type real_type; 178 typedef typename type_traits<element_type>::precision_type precision_type; 179 180 static const unsigned plus_complexity = type_traits<element_type>::plus_complexity; 181 static const unsigned multiplies_complexity = type_traits<element_type>::multiplies_complexity; 182 183 static 184 BOOST_UBLAS_INLINE realboost::numeric::ublas::type_traits185 real_type real (const_reference t) { 186 return type_traits<element_type>::real (t); 187 } 188 static 189 BOOST_UBLAS_INLINE imagboost::numeric::ublas::type_traits190 real_type imag (const_reference t) { 191 return type_traits<element_type>::imag (t); 192 } 193 static 194 BOOST_UBLAS_INLINE conjboost::numeric::ublas::type_traits195 value_type conj (const_reference t) { 196 return type_traits<element_type>::conj (t); 197 } 198 199 static 200 BOOST_UBLAS_INLINE type_absboost::numeric::ublas::type_traits201 real_type type_abs (const_reference t) { 202 return type_traits<element_type>::type_abs (t); 203 } 204 static 205 BOOST_UBLAS_INLINE type_sqrtboost::numeric::ublas::type_traits206 value_type type_sqrt (const_reference t) { 207 return type_traits<element_type>::type_sqrt (t); 208 } 209 210 static 211 BOOST_UBLAS_INLINE norm_1boost::numeric::ublas::type_traits212 real_type norm_1 (const_reference t) { 213 return type_traits<element_type>::norm_1 (t); 214 } 215 static 216 BOOST_UBLAS_INLINE norm_2boost::numeric::ublas::type_traits217 real_type norm_2 (const_reference t) { 218 return type_traits<element_type>::norm_2 (t); 219 } 220 static 221 BOOST_UBLAS_INLINE norm_infboost::numeric::ublas::type_traits222 real_type norm_inf (const_reference t) { 223 return type_traits<element_type>::norm_inf (t); 224 } 225 226 static 227 BOOST_UBLAS_INLINE equalsboost::numeric::ublas::type_traits228 bool equals (const_reference t1, const_reference t2) { 229 return type_traits<element_type>::equals (t1, t2); 230 } 231 }; 232 233 template<class M1, class T2> 234 struct promote_traits<sparse_matrix_element<M1>, T2> { 235 typedef typename promote_traits<typename sparse_matrix_element<M1>::value_type, T2>::promote_type promote_type; 236 }; 237 template<class T1, class M2> 238 struct promote_traits<T1, sparse_matrix_element<M2> > { 239 typedef typename promote_traits<T1, typename sparse_matrix_element<M2>::value_type>::promote_type promote_type; 240 }; 241 template<class M1, class M2> 242 struct promote_traits<sparse_matrix_element<M1>, sparse_matrix_element<M2> > { 243 typedef typename promote_traits<typename sparse_matrix_element<M1>::value_type, 244 typename sparse_matrix_element<M2>::value_type>::promote_type promote_type; 245 }; 246 247 #endif 248 249 /** \brief Index map based sparse matrix of values of type \c T 250 * 251 * This class represents a matrix by using a \c key to value mapping. The default type is 252 * \code template<class T, class L = row_major, class A = map_std<std::size_t, T> > class mapped_matrix; \endcode 253 * So, by default a STL map container is used to associate keys and values. The key is computed depending on 254 * the layout type \c L as \code key = layout_type::element(i, size1_, j, size2_); \endcode 255 * which means \code key = (i*size2+j) \endcode for a row major matrix. 256 * Limitations: The matrix size must not exceed \f$(size1*size2) < \f$ \code std::limits<std::size_t> \endcode. 257 * The \ref find1() and \ref find2() operations have a complexity of at least \f$\mathcal{O}(log(nnz))\f$, depending 258 * on the efficiency of \c std::lower_bound on the key set of the map. 259 * Orientation and storage can also be specified, otherwise a row major orientation is used. 260 * It is \b not required by the storage to initialize elements of the matrix. By default, the orientation is \c row_major. 261 * 262 * \sa fwd.hpp, storage_sparse.hpp 263 * 264 * \tparam T the type of object stored in the matrix (like double, float, complex, etc...) 265 * \tparam L the storage organization. It can be either \c row_major or \c column_major. By default it is \c row_major 266 */ 267 template<class T, class L, class A> 268 class mapped_matrix: 269 public matrix_container<mapped_matrix<T, L, A> > { 270 271 typedef T &true_reference; 272 typedef T *pointer; 273 typedef const T * const_pointer; 274 typedef L layout_type; 275 typedef mapped_matrix<T, L, A> self_type; 276 public: 277 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS 278 using matrix_container<self_type>::operator (); 279 #endif 280 typedef typename A::size_type size_type; 281 typedef typename A::difference_type difference_type; 282 typedef T value_type; 283 typedef A array_type; 284 typedef const T &const_reference; 285 #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE 286 typedef typename detail::map_traits<A, T>::reference reference; 287 #else 288 typedef sparse_matrix_element<self_type> reference; 289 #endif 290 typedef const matrix_reference<const self_type> const_closure_type; 291 typedef matrix_reference<self_type> closure_type; 292 typedef mapped_vector<T, A> vector_temporary_type; 293 typedef self_type matrix_temporary_type; 294 typedef sparse_tag storage_category; 295 typedef typename L::orientation_category orientation_category; 296 297 // Construction and destruction 298 BOOST_UBLAS_INLINE mapped_matrix()299 mapped_matrix (): 300 matrix_container<self_type> (), 301 size1_ (0), size2_ (0), data_ () {} 302 BOOST_UBLAS_INLINE mapped_matrix(size_type size1,size_type size2,size_type non_zeros=0)303 mapped_matrix (size_type size1, size_type size2, size_type non_zeros = 0): 304 matrix_container<self_type> (), 305 size1_ (size1), size2_ (size2), data_ () { 306 detail::map_reserve (data (), restrict_capacity (non_zeros)); 307 } 308 BOOST_UBLAS_INLINE mapped_matrix(const mapped_matrix & m)309 mapped_matrix (const mapped_matrix &m): 310 matrix_container<self_type> (), 311 size1_ (m.size1_), size2_ (m.size2_), data_ (m.data_) {} 312 template<class AE> 313 BOOST_UBLAS_INLINE mapped_matrix(const matrix_expression<AE> & ae,size_type non_zeros=0)314 mapped_matrix (const matrix_expression<AE> &ae, size_type non_zeros = 0): 315 matrix_container<self_type> (), 316 size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), data_ () { 317 detail::map_reserve (data (), restrict_capacity (non_zeros)); 318 matrix_assign<scalar_assign> (*this, ae); 319 } 320 321 // Accessors 322 BOOST_UBLAS_INLINE size1() const323 size_type size1 () const { 324 return size1_; 325 } 326 BOOST_UBLAS_INLINE size2() const327 size_type size2 () const { 328 return size2_; 329 } 330 BOOST_UBLAS_INLINE nnz_capacity() const331 size_type nnz_capacity () const { 332 return detail::map_capacity (data ()); 333 } 334 BOOST_UBLAS_INLINE nnz() const335 size_type nnz () const { 336 return data (). size (); 337 } 338 339 // Storage accessors 340 BOOST_UBLAS_INLINE data() const341 const array_type &data () const { 342 return data_; 343 } 344 BOOST_UBLAS_INLINE data()345 array_type &data () { 346 return data_; 347 } 348 349 // Resizing 350 private: 351 BOOST_UBLAS_INLINE restrict_capacity(size_type non_zeros) const352 size_type restrict_capacity (size_type non_zeros) const { 353 // Guarding against overflow - thanks to Alexei Novakov for the hint. 354 // non_zeros = (std::min) (non_zeros, size1_ * size2_); 355 if (size1_ > 0 && non_zeros / size1_ >= size2_) 356 non_zeros = size1_ * size2_; 357 return non_zeros; 358 } 359 public: 360 BOOST_UBLAS_INLINE resize(size_type size1,size_type size2,bool preserve=true)361 void resize (size_type size1, size_type size2, bool preserve = true) { 362 // FIXME preserve unimplemented 363 BOOST_UBLAS_CHECK (!preserve, internal_logic ()); 364 size1_ = size1; 365 size2_ = size2; 366 data ().clear (); 367 } 368 369 // Reserving 370 BOOST_UBLAS_INLINE reserve(size_type non_zeros,bool preserve=true)371 void reserve (size_type non_zeros, bool preserve = true) { 372 detail::map_reserve (data (), restrict_capacity (non_zeros)); 373 } 374 375 // Element support 376 BOOST_UBLAS_INLINE find_element(size_type i,size_type j)377 pointer find_element (size_type i, size_type j) { 378 return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i, j)); 379 } 380 BOOST_UBLAS_INLINE find_element(size_type i,size_type j) const381 const_pointer find_element (size_type i, size_type j) const { 382 const size_type element = layout_type::element (i, size1_, j, size2_); 383 const_subiterator_type it (data ().find (element)); 384 if (it == data ().end ()) 385 return 0; 386 BOOST_UBLAS_CHECK ((*it).first == element, internal_logic ()); // broken map 387 return &(*it).second; 388 } 389 390 // Element access 391 BOOST_UBLAS_INLINE operator ()(size_type i,size_type j) const392 const_reference operator () (size_type i, size_type j) const { 393 const size_type element = layout_type::element (i, size1_, j, size2_); 394 const_subiterator_type it (data ().find (element)); 395 if (it == data ().end ()) 396 return zero_; 397 BOOST_UBLAS_CHECK ((*it).first == element, internal_logic ()); // broken map 398 return (*it).second; 399 } 400 BOOST_UBLAS_INLINE operator ()(size_type i,size_type j)401 reference operator () (size_type i, size_type j) { 402 #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE 403 const size_type element = layout_type::element (i, size1_, j, size2_); 404 std::pair<subiterator_type, bool> ii (data ().insert (typename array_type::value_type (element, value_type/*zero*/()))); 405 BOOST_UBLAS_CHECK ((ii.first)->first == element, internal_logic ()); // broken map 406 return (ii.first)->second; 407 #else 408 return reference (*this, i, j); 409 #endif 410 } 411 412 // Element assingment 413 BOOST_UBLAS_INLINE insert_element(size_type i,size_type j,const_reference t)414 true_reference insert_element (size_type i, size_type j, const_reference t) { 415 BOOST_UBLAS_CHECK (!find_element (i, j), bad_index ()); // duplicate element 416 const size_type element = layout_type::element (i, size1_, j, size2_); 417 std::pair<subiterator_type, bool> ii (data ().insert (typename array_type::value_type (element, t))); 418 BOOST_UBLAS_CHECK ((ii.first)->first == element, internal_logic ()); // broken map 419 if (!ii.second) // existing element 420 (ii.first)->second = t; 421 return (ii.first)->second; 422 } 423 BOOST_UBLAS_INLINE erase_element(size_type i,size_type j)424 void erase_element (size_type i, size_type j) { 425 subiterator_type it = data ().find (layout_type::element (i, size1_, j, size2_)); 426 if (it == data ().end ()) 427 return; 428 data ().erase (it); 429 } 430 431 // Zeroing 432 BOOST_UBLAS_INLINE clear()433 void clear () { 434 data ().clear (); 435 } 436 437 // Assignment 438 BOOST_UBLAS_INLINE operator =(const mapped_matrix & m)439 mapped_matrix &operator = (const mapped_matrix &m) { 440 if (this != &m) { 441 size1_ = m.size1_; 442 size2_ = m.size2_; 443 data () = m.data (); 444 } 445 return *this; 446 } 447 template<class C> // Container assignment without temporary 448 BOOST_UBLAS_INLINE operator =(const matrix_container<C> & m)449 mapped_matrix &operator = (const matrix_container<C> &m) { 450 resize (m ().size1 (), m ().size2 (), false); 451 assign (m); 452 return *this; 453 } 454 BOOST_UBLAS_INLINE assign_temporary(mapped_matrix & m)455 mapped_matrix &assign_temporary (mapped_matrix &m) { 456 swap (m); 457 return *this; 458 } 459 template<class AE> 460 BOOST_UBLAS_INLINE operator =(const matrix_expression<AE> & ae)461 mapped_matrix &operator = (const matrix_expression<AE> &ae) { 462 self_type temporary (ae, detail::map_capacity (data ())); 463 return assign_temporary (temporary); 464 } 465 template<class AE> 466 BOOST_UBLAS_INLINE assign(const matrix_expression<AE> & ae)467 mapped_matrix &assign (const matrix_expression<AE> &ae) { 468 matrix_assign<scalar_assign> (*this, ae); 469 return *this; 470 } 471 template<class AE> 472 BOOST_UBLAS_INLINE operator +=(const matrix_expression<AE> & ae)473 mapped_matrix& operator += (const matrix_expression<AE> &ae) { 474 self_type temporary (*this + ae, detail::map_capacity (data ())); 475 return assign_temporary (temporary); 476 } 477 template<class C> // Container assignment without temporary 478 BOOST_UBLAS_INLINE operator +=(const matrix_container<C> & m)479 mapped_matrix &operator += (const matrix_container<C> &m) { 480 plus_assign (m); 481 return *this; 482 } 483 template<class AE> 484 BOOST_UBLAS_INLINE plus_assign(const matrix_expression<AE> & ae)485 mapped_matrix &plus_assign (const matrix_expression<AE> &ae) { 486 matrix_assign<scalar_plus_assign> (*this, ae); 487 return *this; 488 } 489 template<class AE> 490 BOOST_UBLAS_INLINE operator -=(const matrix_expression<AE> & ae)491 mapped_matrix& operator -= (const matrix_expression<AE> &ae) { 492 self_type temporary (*this - ae, detail::map_capacity (data ())); 493 return assign_temporary (temporary); 494 } 495 template<class C> // Container assignment without temporary 496 BOOST_UBLAS_INLINE operator -=(const matrix_container<C> & m)497 mapped_matrix &operator -= (const matrix_container<C> &m) { 498 minus_assign (m); 499 return *this; 500 } 501 template<class AE> 502 BOOST_UBLAS_INLINE minus_assign(const matrix_expression<AE> & ae)503 mapped_matrix &minus_assign (const matrix_expression<AE> &ae) { 504 matrix_assign<scalar_minus_assign> (*this, ae); 505 return *this; 506 } 507 template<class AT> 508 BOOST_UBLAS_INLINE operator *=(const AT & at)509 mapped_matrix& operator *= (const AT &at) { 510 matrix_assign_scalar<scalar_multiplies_assign> (*this, at); 511 return *this; 512 } 513 template<class AT> 514 BOOST_UBLAS_INLINE operator /=(const AT & at)515 mapped_matrix& operator /= (const AT &at) { 516 matrix_assign_scalar<scalar_divides_assign> (*this, at); 517 return *this; 518 } 519 520 // Swapping 521 BOOST_UBLAS_INLINE swap(mapped_matrix & m)522 void swap (mapped_matrix &m) { 523 if (this != &m) { 524 std::swap (size1_, m.size1_); 525 std::swap (size2_, m.size2_); 526 data ().swap (m.data ()); 527 } 528 } 529 BOOST_UBLAS_INLINE swap(mapped_matrix & m1,mapped_matrix & m2)530 friend void swap (mapped_matrix &m1, mapped_matrix &m2) { 531 m1.swap (m2); 532 } 533 534 // Iterator types 535 private: 536 // Use storage iterator 537 typedef typename A::const_iterator const_subiterator_type; 538 typedef typename A::iterator subiterator_type; 539 540 BOOST_UBLAS_INLINE at_element(size_type i,size_type j)541 true_reference at_element (size_type i, size_type j) { 542 const size_type element = layout_type::element (i, size1_, j, size2_); 543 subiterator_type it (data ().find (element)); 544 BOOST_UBLAS_CHECK (it != data ().end(), bad_index ()); 545 BOOST_UBLAS_CHECK ((*it).first == element, internal_logic ()); // broken map 546 return it->second; 547 } 548 549 public: 550 class const_iterator1; 551 class iterator1; 552 class const_iterator2; 553 class iterator2; 554 typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1; 555 typedef reverse_iterator_base1<iterator1> reverse_iterator1; 556 typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2; 557 typedef reverse_iterator_base2<iterator2> reverse_iterator2; 558 559 // Element lookup 560 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. find1(int rank,size_type i,size_type j,int direction=1) const561 const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const { 562 const_subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_))); 563 const_subiterator_type it_end (data ().end ()); 564 size_type index1 = size_type (-1); 565 size_type index2 = size_type (-1); 566 while (rank == 1 && it != it_end) { 567 index1 = layout_type::index_i ((*it).first, size1_, size2_); 568 index2 = layout_type::index_j ((*it).first, size1_, size2_); 569 if (direction > 0) { 570 if ((index1 >= i && index2 == j) || (i >= size1_)) 571 break; 572 ++ i; 573 } else /* if (direction < 0) */ { 574 if ((index1 <= i && index2 == j) || (i == 0)) 575 break; 576 -- i; 577 } 578 it = data ().lower_bound (layout_type::address (i, size1_, j, size2_)); 579 } 580 if (rank == 1 && index2 != j) { 581 if (direction > 0) 582 i = size1_; 583 else /* if (direction < 0) */ 584 i = 0; 585 rank = 0; 586 } 587 return const_iterator1 (*this, rank, i, j, it); 588 } 589 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. find1(int rank,size_type i,size_type j,int direction=1)590 iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) { 591 subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_))); 592 subiterator_type it_end (data ().end ()); 593 size_type index1 = size_type (-1); 594 size_type index2 = size_type (-1); 595 while (rank == 1 && it != it_end) { 596 index1 = layout_type::index_i ((*it).first, size1_, size2_); 597 index2 = layout_type::index_j ((*it).first, size1_, size2_); 598 if (direction > 0) { 599 if ((index1 >= i && index2 == j) || (i >= size1_)) 600 break; 601 ++ i; 602 } else /* if (direction < 0) */ { 603 if ((index1 <= i && index2 == j) || (i == 0)) 604 break; 605 -- i; 606 } 607 it = data ().lower_bound (layout_type::address (i, size1_, j, size2_)); 608 } 609 if (rank == 1 && index2 != j) { 610 if (direction > 0) 611 i = size1_; 612 else /* if (direction < 0) */ 613 i = 0; 614 rank = 0; 615 } 616 return iterator1 (*this, rank, i, j, it); 617 } 618 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. find2(int rank,size_type i,size_type j,int direction=1) const619 const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const { 620 const_subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_))); 621 const_subiterator_type it_end (data ().end ()); 622 size_type index1 = size_type (-1); 623 size_type index2 = size_type (-1); 624 while (rank == 1 && it != it_end) { 625 index1 = layout_type::index_i ((*it).first, size1_, size2_); 626 index2 = layout_type::index_j ((*it).first, size1_, size2_); 627 if (direction > 0) { 628 if ((index2 >= j && index1 == i) || (j >= size2_)) 629 break; 630 ++ j; 631 } else /* if (direction < 0) */ { 632 if ((index2 <= j && index1 == i) || (j == 0)) 633 break; 634 -- j; 635 } 636 it = data ().lower_bound (layout_type::address (i, size1_, j, size2_)); 637 } 638 if (rank == 1 && index1 != i) { 639 if (direction > 0) 640 j = size2_; 641 else /* if (direction < 0) */ 642 j = 0; 643 rank = 0; 644 } 645 return const_iterator2 (*this, rank, i, j, it); 646 } 647 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. find2(int rank,size_type i,size_type j,int direction=1)648 iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) { 649 subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_))); 650 subiterator_type it_end (data ().end ()); 651 size_type index1 = size_type (-1); 652 size_type index2 = size_type (-1); 653 while (rank == 1 && it != it_end) { 654 index1 = layout_type::index_i ((*it).first, size1_, size2_); 655 index2 = layout_type::index_j ((*it).first, size1_, size2_); 656 if (direction > 0) { 657 if ((index2 >= j && index1 == i) || (j >= size2_)) 658 break; 659 ++ j; 660 } else /* if (direction < 0) */ { 661 if ((index2 <= j && index1 == i) || (j == 0)) 662 break; 663 -- j; 664 } 665 it = data ().lower_bound (layout_type::address (i, size1_, j, size2_)); 666 } 667 if (rank == 1 && index1 != i) { 668 if (direction > 0) 669 j = size2_; 670 else /* if (direction < 0) */ 671 j = 0; 672 rank = 0; 673 } 674 return iterator2 (*this, rank, i, j, it); 675 } 676 677 678 class const_iterator1: 679 public container_const_reference<mapped_matrix>, 680 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag, 681 const_iterator1, value_type> { 682 public: 683 typedef typename mapped_matrix::value_type value_type; 684 typedef typename mapped_matrix::difference_type difference_type; 685 typedef typename mapped_matrix::const_reference reference; 686 typedef const typename mapped_matrix::pointer pointer; 687 688 typedef const_iterator2 dual_iterator_type; 689 typedef const_reverse_iterator2 dual_reverse_iterator_type; 690 691 // Construction and destruction 692 BOOST_UBLAS_INLINE const_iterator1()693 const_iterator1 (): 694 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), it_ () {} 695 BOOST_UBLAS_INLINE const_iterator1(const self_type & m,int rank,size_type i,size_type j,const const_subiterator_type & it)696 const_iterator1 (const self_type &m, int rank, size_type i, size_type j, const const_subiterator_type &it): 697 container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), it_ (it) {} 698 BOOST_UBLAS_INLINE const_iterator1(const iterator1 & it)699 const_iterator1 (const iterator1 &it): 700 container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), it_ (it.it_) {} 701 702 // Arithmetic 703 BOOST_UBLAS_INLINE operator ++()704 const_iterator1 &operator ++ () { 705 if (rank_ == 1 && layout_type::fast_i ()) 706 ++ it_; 707 else 708 *this = (*this) ().find1 (rank_, index1 () + 1, j_, 1); 709 return *this; 710 } 711 BOOST_UBLAS_INLINE operator --()712 const_iterator1 &operator -- () { 713 if (rank_ == 1 && layout_type::fast_i ()) 714 -- it_; 715 else 716 *this = (*this) ().find1 (rank_, index1 () - 1, j_, -1); 717 return *this; 718 } 719 720 // Dereference 721 BOOST_UBLAS_INLINE operator *() const722 const_reference operator * () const { 723 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ()); 724 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ()); 725 if (rank_ == 1) { 726 return (*it_).second; 727 } else { 728 return (*this) () (i_, j_); 729 } 730 } 731 732 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 733 BOOST_UBLAS_INLINE 734 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 735 typename self_type:: 736 #endif begin() const737 const_iterator2 begin () const { 738 const self_type &m = (*this) (); 739 return m.find2 (1, index1 (), 0); 740 } 741 BOOST_UBLAS_INLINE 742 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 743 typename self_type:: 744 #endif cbegin() const745 const_iterator2 cbegin () const { 746 return begin (); 747 } 748 BOOST_UBLAS_INLINE 749 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 750 typename self_type:: 751 #endif end() const752 const_iterator2 end () const { 753 const self_type &m = (*this) (); 754 return m.find2 (1, index1 (), m.size2 ()); 755 } 756 BOOST_UBLAS_INLINE 757 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 758 typename self_type:: 759 #endif cend() const760 const_iterator2 cend () const { 761 return end (); 762 } 763 BOOST_UBLAS_INLINE 764 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 765 typename self_type:: 766 #endif rbegin() const767 const_reverse_iterator2 rbegin () const { 768 return const_reverse_iterator2 (end ()); 769 } 770 BOOST_UBLAS_INLINE 771 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 772 typename self_type:: 773 #endif crbegin() const774 const_reverse_iterator2 crbegin () const { 775 return rbegin (); 776 } 777 BOOST_UBLAS_INLINE 778 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 779 typename self_type:: 780 #endif rend() const781 const_reverse_iterator2 rend () const { 782 return const_reverse_iterator2 (begin ()); 783 } 784 BOOST_UBLAS_INLINE 785 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 786 typename self_type:: 787 #endif crend() const788 const_reverse_iterator2 crend () const { 789 return rend (); 790 } 791 #endif 792 793 // Indices 794 BOOST_UBLAS_INLINE index1() const795 size_type index1 () const { 796 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ()); 797 if (rank_ == 1) { 798 const self_type &m = (*this) (); 799 BOOST_UBLAS_CHECK (layout_type::index_i ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size1 (), bad_index ()); 800 return layout_type::index_i ((*it_).first, m.size1 (), m.size2 ()); 801 } else { 802 return i_; 803 } 804 } 805 BOOST_UBLAS_INLINE index2() const806 size_type index2 () const { 807 if (rank_ == 1) { 808 const self_type &m = (*this) (); 809 BOOST_UBLAS_CHECK (layout_type::index_j ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size2 (), bad_index ()); 810 return layout_type::index_j ((*it_).first, m.size1 (), m.size2 ()); 811 } else { 812 return j_; 813 } 814 } 815 816 // Assignment 817 BOOST_UBLAS_INLINE operator =(const const_iterator1 & it)818 const_iterator1 &operator = (const const_iterator1 &it) { 819 container_const_reference<self_type>::assign (&it ()); 820 rank_ = it.rank_; 821 i_ = it.i_; 822 j_ = it.j_; 823 it_ = it.it_; 824 return *this; 825 } 826 827 // Comparison 828 BOOST_UBLAS_INLINE operator ==(const const_iterator1 & it) const829 bool operator == (const const_iterator1 &it) const { 830 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 831 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ()); 832 if (rank_ == 1 || it.rank_ == 1) { 833 return it_ == it.it_; 834 } else { 835 return i_ == it.i_ && j_ == it.j_; 836 } 837 } 838 839 private: 840 int rank_; 841 size_type i_; 842 size_type j_; 843 const_subiterator_type it_; 844 }; 845 846 BOOST_UBLAS_INLINE begin1() const847 const_iterator1 begin1 () const { 848 return find1 (0, 0, 0); 849 } 850 BOOST_UBLAS_INLINE cbegin1() const851 const_iterator1 cbegin1 () const { 852 return begin1 (); 853 } 854 BOOST_UBLAS_INLINE end1() const855 const_iterator1 end1 () const { 856 return find1 (0, size1_, 0); 857 } 858 BOOST_UBLAS_INLINE cend1() const859 const_iterator1 cend1 () const { 860 return end1 (); 861 } 862 863 class iterator1: 864 public container_reference<mapped_matrix>, 865 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag, 866 iterator1, value_type> { 867 public: 868 typedef typename mapped_matrix::value_type value_type; 869 typedef typename mapped_matrix::difference_type difference_type; 870 typedef typename mapped_matrix::true_reference reference; 871 typedef typename mapped_matrix::pointer pointer; 872 873 typedef iterator2 dual_iterator_type; 874 typedef reverse_iterator2 dual_reverse_iterator_type; 875 876 // Construction and destruction 877 BOOST_UBLAS_INLINE iterator1()878 iterator1 (): 879 container_reference<self_type> (), rank_ (), i_ (), j_ (), it_ () {} 880 BOOST_UBLAS_INLINE iterator1(self_type & m,int rank,size_type i,size_type j,const subiterator_type & it)881 iterator1 (self_type &m, int rank, size_type i, size_type j, const subiterator_type &it): 882 container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), it_ (it) {} 883 884 // Arithmetic 885 BOOST_UBLAS_INLINE operator ++()886 iterator1 &operator ++ () { 887 if (rank_ == 1 && layout_type::fast_i ()) 888 ++ it_; 889 else 890 *this = (*this) ().find1 (rank_, index1 () + 1, j_, 1); 891 return *this; 892 } 893 BOOST_UBLAS_INLINE operator --()894 iterator1 &operator -- () { 895 if (rank_ == 1 && layout_type::fast_i ()) 896 -- it_; 897 else 898 *this = (*this) ().find1 (rank_, index1 () - 1, j_, -1); 899 return *this; 900 } 901 902 // Dereference 903 BOOST_UBLAS_INLINE operator *() const904 reference operator * () const { 905 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ()); 906 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ()); 907 if (rank_ == 1) { 908 return (*it_).second; 909 } else { 910 return (*this) ().at_element (i_, j_); 911 } 912 } 913 914 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 915 BOOST_UBLAS_INLINE 916 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 917 typename self_type:: 918 #endif begin() const919 iterator2 begin () const { 920 self_type &m = (*this) (); 921 return m.find2 (1, index1 (), 0); 922 } 923 BOOST_UBLAS_INLINE 924 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 925 typename self_type:: 926 #endif end() const927 iterator2 end () const { 928 self_type &m = (*this) (); 929 return m.find2 (1, index1 (), m.size2 ()); 930 } 931 BOOST_UBLAS_INLINE 932 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 933 typename self_type:: 934 #endif rbegin() const935 reverse_iterator2 rbegin () const { 936 return reverse_iterator2 (end ()); 937 } 938 BOOST_UBLAS_INLINE 939 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 940 typename self_type:: 941 #endif rend() const942 reverse_iterator2 rend () const { 943 return reverse_iterator2 (begin ()); 944 } 945 #endif 946 947 // Indices 948 BOOST_UBLAS_INLINE index1() const949 size_type index1 () const { 950 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ()); 951 if (rank_ == 1) { 952 const self_type &m = (*this) (); 953 BOOST_UBLAS_CHECK (layout_type::index_i ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size1 (), bad_index ()); 954 return layout_type::index_i ((*it_).first, m.size1 (), m.size2 ()); 955 } else { 956 return i_; 957 } 958 } 959 BOOST_UBLAS_INLINE index2() const960 size_type index2 () const { 961 if (rank_ == 1) { 962 const self_type &m = (*this) (); 963 BOOST_UBLAS_CHECK (layout_type::index_j ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size2 (), bad_index ()); 964 return layout_type::index_j ((*it_).first, m.size1 (), m.size2 ()); 965 } else { 966 return j_; 967 } 968 } 969 970 // Assignment 971 BOOST_UBLAS_INLINE operator =(const iterator1 & it)972 iterator1 &operator = (const iterator1 &it) { 973 container_reference<self_type>::assign (&it ()); 974 rank_ = it.rank_; 975 i_ = it.i_; 976 j_ = it.j_; 977 it_ = it.it_; 978 return *this; 979 } 980 981 // Comparison 982 BOOST_UBLAS_INLINE operator ==(const iterator1 & it) const983 bool operator == (const iterator1 &it) const { 984 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 985 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ()); 986 if (rank_ == 1 || it.rank_ == 1) { 987 return it_ == it.it_; 988 } else { 989 return i_ == it.i_ && j_ == it.j_; 990 } 991 } 992 993 private: 994 int rank_; 995 size_type i_; 996 size_type j_; 997 subiterator_type it_; 998 999 friend class const_iterator1; 1000 }; 1001 1002 BOOST_UBLAS_INLINE begin1()1003 iterator1 begin1 () { 1004 return find1 (0, 0, 0); 1005 } 1006 BOOST_UBLAS_INLINE end1()1007 iterator1 end1 () { 1008 return find1 (0, size1_, 0); 1009 } 1010 1011 class const_iterator2: 1012 public container_const_reference<mapped_matrix>, 1013 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag, 1014 const_iterator2, value_type> { 1015 public: 1016 typedef typename mapped_matrix::value_type value_type; 1017 typedef typename mapped_matrix::difference_type difference_type; 1018 typedef typename mapped_matrix::const_reference reference; 1019 typedef const typename mapped_matrix::pointer pointer; 1020 1021 typedef const_iterator1 dual_iterator_type; 1022 typedef const_reverse_iterator1 dual_reverse_iterator_type; 1023 1024 // Construction and destruction 1025 BOOST_UBLAS_INLINE const_iterator2()1026 const_iterator2 (): 1027 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), it_ () {} 1028 BOOST_UBLAS_INLINE const_iterator2(const self_type & m,int rank,size_type i,size_type j,const const_subiterator_type & it)1029 const_iterator2 (const self_type &m, int rank, size_type i, size_type j, const const_subiterator_type &it): 1030 container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), it_ (it) {} 1031 BOOST_UBLAS_INLINE const_iterator2(const iterator2 & it)1032 const_iterator2 (const iterator2 &it): 1033 container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), it_ (it.it_) {} 1034 1035 // Arithmetic 1036 BOOST_UBLAS_INLINE operator ++()1037 const_iterator2 &operator ++ () { 1038 if (rank_ == 1 && layout_type::fast_j ()) 1039 ++ it_; 1040 else 1041 *this = (*this) ().find2 (rank_, i_, index2 () + 1, 1); 1042 return *this; 1043 } 1044 BOOST_UBLAS_INLINE operator --()1045 const_iterator2 &operator -- () { 1046 if (rank_ == 1 && layout_type::fast_j ()) 1047 -- it_; 1048 else 1049 *this = (*this) ().find2 (rank_, i_, index2 () - 1, -1); 1050 return *this; 1051 } 1052 1053 // Dereference 1054 BOOST_UBLAS_INLINE operator *() const1055 const_reference operator * () const { 1056 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ()); 1057 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ()); 1058 if (rank_ == 1) { 1059 return (*it_).second; 1060 } else { 1061 return (*this) () (i_, j_); 1062 } 1063 } 1064 1065 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 1066 BOOST_UBLAS_INLINE 1067 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 1068 typename self_type:: 1069 #endif begin() const1070 const_iterator1 begin () const { 1071 const self_type &m = (*this) (); 1072 return m.find1 (1, 0, index2 ()); 1073 } 1074 BOOST_UBLAS_INLINE 1075 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 1076 typename self_type:: 1077 #endif cbegin() const1078 const_iterator1 cbegin () const { 1079 return begin (); 1080 } 1081 BOOST_UBLAS_INLINE 1082 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 1083 typename self_type:: 1084 #endif end() const1085 const_iterator1 end () const { 1086 const self_type &m = (*this) (); 1087 return m.find1 (1, m.size1 (), index2 ()); 1088 } 1089 BOOST_UBLAS_INLINE 1090 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 1091 typename self_type:: 1092 #endif cend() const1093 const_iterator1 cend () const { 1094 return end (); 1095 } 1096 BOOST_UBLAS_INLINE 1097 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 1098 typename self_type:: 1099 #endif rbegin() const1100 const_reverse_iterator1 rbegin () const { 1101 return const_reverse_iterator1 (end ()); 1102 } 1103 BOOST_UBLAS_INLINE 1104 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 1105 typename self_type:: 1106 #endif crbegin() const1107 const_reverse_iterator1 crbegin () const { 1108 return rbegin (); 1109 } 1110 BOOST_UBLAS_INLINE 1111 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 1112 typename self_type:: 1113 #endif rend() const1114 const_reverse_iterator1 rend () const { 1115 return const_reverse_iterator1 (begin ()); 1116 } 1117 BOOST_UBLAS_INLINE 1118 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 1119 typename self_type:: 1120 #endif crend() const1121 const_reverse_iterator1 crend () const { 1122 return rend (); 1123 } 1124 #endif 1125 1126 // Indices 1127 BOOST_UBLAS_INLINE index1() const1128 size_type index1 () const { 1129 if (rank_ == 1) { 1130 const self_type &m = (*this) (); 1131 BOOST_UBLAS_CHECK (layout_type::index_i ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size1 (), bad_index ()); 1132 return layout_type::index_i ((*it_).first, m.size1 (), m.size2 ()); 1133 } else { 1134 return i_; 1135 } 1136 } 1137 BOOST_UBLAS_INLINE index2() const1138 size_type index2 () const { 1139 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ()); 1140 if (rank_ == 1) { 1141 const self_type &m = (*this) (); 1142 BOOST_UBLAS_CHECK (layout_type::index_j ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size2 (), bad_index ()); 1143 return layout_type::index_j ((*it_).first, m.size1 (), m.size2 ()); 1144 } else { 1145 return j_; 1146 } 1147 } 1148 1149 // Assignment 1150 BOOST_UBLAS_INLINE operator =(const const_iterator2 & it)1151 const_iterator2 &operator = (const const_iterator2 &it) { 1152 container_const_reference<self_type>::assign (&it ()); 1153 rank_ = it.rank_; 1154 i_ = it.i_; 1155 j_ = it.j_; 1156 it_ = it.it_; 1157 return *this; 1158 } 1159 1160 // Comparison 1161 BOOST_UBLAS_INLINE operator ==(const const_iterator2 & it) const1162 bool operator == (const const_iterator2 &it) const { 1163 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 1164 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ()); 1165 if (rank_ == 1 || it.rank_ == 1) { 1166 return it_ == it.it_; 1167 } else { 1168 return i_ == it.i_ && j_ == it.j_; 1169 } 1170 } 1171 1172 private: 1173 int rank_; 1174 size_type i_; 1175 size_type j_; 1176 const_subiterator_type it_; 1177 }; 1178 1179 BOOST_UBLAS_INLINE begin2() const1180 const_iterator2 begin2 () const { 1181 return find2 (0, 0, 0); 1182 } 1183 BOOST_UBLAS_INLINE cbegin2() const1184 const_iterator2 cbegin2 () const { 1185 return begin2 (); 1186 } 1187 BOOST_UBLAS_INLINE end2() const1188 const_iterator2 end2 () const { 1189 return find2 (0, 0, size2_); 1190 } 1191 BOOST_UBLAS_INLINE cend2() const1192 const_iterator2 cend2 () const { 1193 return end2 (); 1194 } 1195 1196 class iterator2: 1197 public container_reference<mapped_matrix>, 1198 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag, 1199 iterator2, value_type> { 1200 public: 1201 typedef typename mapped_matrix::value_type value_type; 1202 typedef typename mapped_matrix::difference_type difference_type; 1203 typedef typename mapped_matrix::true_reference reference; 1204 typedef typename mapped_matrix::pointer pointer; 1205 1206 typedef iterator1 dual_iterator_type; 1207 typedef reverse_iterator1 dual_reverse_iterator_type; 1208 1209 // Construction and destruction 1210 BOOST_UBLAS_INLINE iterator2()1211 iterator2 (): 1212 container_reference<self_type> (), rank_ (), i_ (), j_ (), it_ () {} 1213 BOOST_UBLAS_INLINE iterator2(self_type & m,int rank,size_type i,size_type j,const subiterator_type & it)1214 iterator2 (self_type &m, int rank, size_type i, size_type j, const subiterator_type &it): 1215 container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), it_ (it) {} 1216 1217 // Arithmetic 1218 BOOST_UBLAS_INLINE operator ++()1219 iterator2 &operator ++ () { 1220 if (rank_ == 1 && layout_type::fast_j ()) 1221 ++ it_; 1222 else 1223 *this = (*this) ().find2 (rank_, i_, index2 () + 1, 1); 1224 return *this; 1225 } 1226 BOOST_UBLAS_INLINE operator --()1227 iterator2 &operator -- () { 1228 if (rank_ == 1 && layout_type::fast_j ()) 1229 -- it_; 1230 else 1231 *this = (*this) ().find2 (rank_, i_, index2 () - 1, -1); 1232 return *this; 1233 } 1234 1235 // Dereference 1236 BOOST_UBLAS_INLINE operator *() const1237 reference operator * () const { 1238 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ()); 1239 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ()); 1240 if (rank_ == 1) { 1241 return (*it_).second; 1242 } else { 1243 return (*this) ().at_element (i_, j_); 1244 } 1245 } 1246 1247 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 1248 BOOST_UBLAS_INLINE 1249 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 1250 typename self_type:: 1251 #endif begin() const1252 iterator1 begin () const { 1253 self_type &m = (*this) (); 1254 return m.find1 (1, 0, index2 ()); 1255 } 1256 BOOST_UBLAS_INLINE 1257 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 1258 typename self_type:: 1259 #endif end() const1260 iterator1 end () const { 1261 self_type &m = (*this) (); 1262 return m.find1 (1, m.size1 (), index2 ()); 1263 } 1264 BOOST_UBLAS_INLINE 1265 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 1266 typename self_type:: 1267 #endif rbegin() const1268 reverse_iterator1 rbegin () const { 1269 return reverse_iterator1 (end ()); 1270 } 1271 BOOST_UBLAS_INLINE 1272 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 1273 typename self_type:: 1274 #endif rend() const1275 reverse_iterator1 rend () const { 1276 return reverse_iterator1 (begin ()); 1277 } 1278 #endif 1279 1280 // Indices 1281 BOOST_UBLAS_INLINE index1() const1282 size_type index1 () const { 1283 if (rank_ == 1) { 1284 const self_type &m = (*this) (); 1285 BOOST_UBLAS_CHECK (layout_type::index_i ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size1 (), bad_index ()); 1286 return layout_type::index_i ((*it_).first, m.size1 (), m.size2 ()); 1287 } else { 1288 return i_; 1289 } 1290 } 1291 BOOST_UBLAS_INLINE index2() const1292 size_type index2 () const { 1293 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ()); 1294 if (rank_ == 1) { 1295 const self_type &m = (*this) (); 1296 BOOST_UBLAS_CHECK (layout_type::index_j ((*it_).first, m.size1 (), m.size2 ()) < (*this) ().size2 (), bad_index ()); 1297 return layout_type::index_j ((*it_).first, m.size1 (), m.size2 ()); 1298 } else { 1299 return j_; 1300 } 1301 } 1302 1303 // Assignment 1304 BOOST_UBLAS_INLINE operator =(const iterator2 & it)1305 iterator2 &operator = (const iterator2 &it) { 1306 container_reference<self_type>::assign (&it ()); 1307 rank_ = it.rank_; 1308 i_ = it.i_; 1309 j_ = it.j_; 1310 it_ = it.it_; 1311 return *this; 1312 } 1313 1314 // Comparison 1315 BOOST_UBLAS_INLINE operator ==(const iterator2 & it) const1316 bool operator == (const iterator2 &it) const { 1317 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 1318 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ()); 1319 if (rank_ == 1 || it.rank_ == 1) { 1320 return it_ == it.it_; 1321 } else { 1322 return i_ == it.i_ && j_ == it.j_; 1323 } 1324 } 1325 1326 private: 1327 int rank_; 1328 size_type i_; 1329 size_type j_; 1330 subiterator_type it_; 1331 1332 friend class const_iterator2; 1333 }; 1334 1335 BOOST_UBLAS_INLINE begin2()1336 iterator2 begin2 () { 1337 return find2 (0, 0, 0); 1338 } 1339 BOOST_UBLAS_INLINE end2()1340 iterator2 end2 () { 1341 return find2 (0, 0, size2_); 1342 } 1343 1344 // Reverse iterators 1345 1346 BOOST_UBLAS_INLINE rbegin1() const1347 const_reverse_iterator1 rbegin1 () const { 1348 return const_reverse_iterator1 (end1 ()); 1349 } 1350 BOOST_UBLAS_INLINE crbegin1() const1351 const_reverse_iterator1 crbegin1 () const { 1352 return rbegin1 (); 1353 } 1354 BOOST_UBLAS_INLINE rend1() const1355 const_reverse_iterator1 rend1 () const { 1356 return const_reverse_iterator1 (begin1 ()); 1357 } 1358 BOOST_UBLAS_INLINE crend1() const1359 const_reverse_iterator1 crend1 () const { 1360 return rend1 (); 1361 } 1362 1363 BOOST_UBLAS_INLINE rbegin1()1364 reverse_iterator1 rbegin1 () { 1365 return reverse_iterator1 (end1 ()); 1366 } 1367 BOOST_UBLAS_INLINE rend1()1368 reverse_iterator1 rend1 () { 1369 return reverse_iterator1 (begin1 ()); 1370 } 1371 1372 BOOST_UBLAS_INLINE rbegin2() const1373 const_reverse_iterator2 rbegin2 () const { 1374 return const_reverse_iterator2 (end2 ()); 1375 } 1376 BOOST_UBLAS_INLINE crbegin2() const1377 const_reverse_iterator2 crbegin2 () const { 1378 return rbegin2 (); 1379 } 1380 BOOST_UBLAS_INLINE rend2() const1381 const_reverse_iterator2 rend2 () const { 1382 return const_reverse_iterator2 (begin2 ()); 1383 } 1384 BOOST_UBLAS_INLINE crend2() const1385 const_reverse_iterator2 crend2 () const { 1386 return rend2 (); 1387 } 1388 1389 BOOST_UBLAS_INLINE rbegin2()1390 reverse_iterator2 rbegin2 () { 1391 return reverse_iterator2 (end2 ()); 1392 } 1393 BOOST_UBLAS_INLINE rend2()1394 reverse_iterator2 rend2 () { 1395 return reverse_iterator2 (begin2 ()); 1396 } 1397 1398 // Serialization 1399 template<class Archive> serialize(Archive & ar,const unsigned int)1400 void serialize(Archive & ar, const unsigned int /* file_version */){ 1401 serialization::collection_size_type s1 (size1_); 1402 serialization::collection_size_type s2 (size2_); 1403 ar & serialization::make_nvp("size1",s1); 1404 ar & serialization::make_nvp("size2",s2); 1405 if (Archive::is_loading::value) { 1406 size1_ = s1; 1407 size2_ = s2; 1408 } 1409 ar & serialization::make_nvp("data", data_); 1410 } 1411 1412 private: 1413 size_type size1_; 1414 size_type size2_; 1415 array_type data_; 1416 static const value_type zero_; 1417 }; 1418 1419 template<class T, class L, class A> 1420 const typename mapped_matrix<T, L, A>::value_type mapped_matrix<T, L, A>::zero_ = value_type/*zero*/(); 1421 1422 1423 // Vector index map based sparse matrix class 1424 template<class T, class L, class A> 1425 class mapped_vector_of_mapped_vector: 1426 public matrix_container<mapped_vector_of_mapped_vector<T, L, A> > { 1427 1428 typedef T &true_reference; 1429 typedef T *pointer; 1430 typedef const T *const_pointer; 1431 typedef A array_type; 1432 typedef const A const_array_type; 1433 typedef L layout_type; 1434 typedef mapped_vector_of_mapped_vector<T, L, A> self_type; 1435 public: 1436 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS 1437 using matrix_container<self_type>::operator (); 1438 #endif 1439 typedef typename A::size_type size_type; 1440 typedef typename A::difference_type difference_type; 1441 typedef T value_type; 1442 typedef const T &const_reference; 1443 #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE 1444 typedef typename detail::map_traits<typename A::data_value_type, T>::reference reference; 1445 #else 1446 typedef sparse_matrix_element<self_type> reference; 1447 #endif 1448 typedef const matrix_reference<const self_type> const_closure_type; 1449 typedef matrix_reference<self_type> closure_type; 1450 typedef mapped_vector<T> vector_temporary_type; 1451 typedef self_type matrix_temporary_type; 1452 typedef typename A::value_type::second_type vector_data_value_type; 1453 typedef sparse_tag storage_category; 1454 typedef typename L::orientation_category orientation_category; 1455 1456 // Construction and destruction 1457 BOOST_UBLAS_INLINE mapped_vector_of_mapped_vector()1458 mapped_vector_of_mapped_vector (): 1459 matrix_container<self_type> (), 1460 size1_ (0), size2_ (0), data_ () { 1461 data_ [layout_type::size_M (size1_, size2_)] = vector_data_value_type (); 1462 } 1463 BOOST_UBLAS_INLINE mapped_vector_of_mapped_vector(size_type size1,size_type size2,size_type non_zeros=0)1464 mapped_vector_of_mapped_vector (size_type size1, size_type size2, size_type non_zeros = 0): 1465 matrix_container<self_type> (), 1466 size1_ (size1), size2_ (size2), data_ () { 1467 data_ [layout_type::size_M (size1_, size2_)] = vector_data_value_type (); 1468 } 1469 BOOST_UBLAS_INLINE mapped_vector_of_mapped_vector(const mapped_vector_of_mapped_vector & m)1470 mapped_vector_of_mapped_vector (const mapped_vector_of_mapped_vector &m): 1471 matrix_container<self_type> (), 1472 size1_ (m.size1_), size2_ (m.size2_), data_ (m.data_) {} 1473 template<class AE> 1474 BOOST_UBLAS_INLINE mapped_vector_of_mapped_vector(const matrix_expression<AE> & ae,size_type non_zeros=0)1475 mapped_vector_of_mapped_vector (const matrix_expression<AE> &ae, size_type non_zeros = 0): 1476 matrix_container<self_type> (), 1477 size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), data_ () { 1478 data_ [layout_type::size_M (size1_, size2_)] = vector_data_value_type (); 1479 matrix_assign<scalar_assign> (*this, ae); 1480 } 1481 1482 // Accessors 1483 BOOST_UBLAS_INLINE size1() const1484 size_type size1 () const { 1485 return size1_; 1486 } 1487 BOOST_UBLAS_INLINE size2() const1488 size_type size2 () const { 1489 return size2_; 1490 } 1491 BOOST_UBLAS_INLINE nnz_capacity() const1492 size_type nnz_capacity () const { 1493 size_type non_zeros = 0; 1494 for (vector_const_subiterator_type itv = data_ ().begin (); itv != data_ ().end (); ++ itv) 1495 non_zeros += detail::map_capacity (*itv); 1496 return non_zeros; 1497 } 1498 BOOST_UBLAS_INLINE nnz() const1499 size_type nnz () const { 1500 size_type filled = 0; 1501 for (vector_const_subiterator_type itv = data_ ().begin (); itv != data_ ().end (); ++ itv) 1502 filled += (*itv).size (); 1503 return filled; 1504 } 1505 1506 // Storage accessors 1507 BOOST_UBLAS_INLINE data() const1508 const_array_type &data () const { 1509 return data_; 1510 } 1511 BOOST_UBLAS_INLINE data()1512 array_type &data () { 1513 return data_; 1514 } 1515 1516 // Resizing 1517 BOOST_UBLAS_INLINE resize(size_type size1,size_type size2,bool preserve=true)1518 void resize (size_type size1, size_type size2, bool preserve = true) { 1519 // FIXME preserve unimplemented 1520 BOOST_UBLAS_CHECK (!preserve, internal_logic ()); 1521 size1_ = size1; 1522 size2_ = size2; 1523 data ().clear (); 1524 data () [layout_type::size_M (size1_, size2_)] = vector_data_value_type (); 1525 } 1526 1527 // Element support 1528 BOOST_UBLAS_INLINE find_element(size_type i,size_type j)1529 pointer find_element (size_type i, size_type j) { 1530 return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i, j)); 1531 } 1532 BOOST_UBLAS_INLINE find_element(size_type i,size_type j) const1533 const_pointer find_element (size_type i, size_type j) const { 1534 const size_type element1 = layout_type::index_M (i, j); 1535 const size_type element2 = layout_type::index_m (i, j); 1536 vector_const_subiterator_type itv (data ().find (element1)); 1537 if (itv == data ().end ()) 1538 return 0; 1539 BOOST_UBLAS_CHECK ((*itv).first == element1, internal_logic ()); // broken map 1540 const_subiterator_type it ((*itv).second.find (element2)); 1541 if (it == (*itv).second.end ()) 1542 return 0; 1543 BOOST_UBLAS_CHECK ((*it).first == element2, internal_logic ()); // broken map 1544 return &(*it).second; 1545 } 1546 1547 // Element access 1548 BOOST_UBLAS_INLINE operator ()(size_type i,size_type j) const1549 const_reference operator () (size_type i, size_type j) const { 1550 const size_type element1 = layout_type::index_M (i, j); 1551 const size_type element2 = layout_type::index_m (i, j); 1552 vector_const_subiterator_type itv (data ().find (element1)); 1553 if (itv == data ().end ()) 1554 return zero_; 1555 BOOST_UBLAS_CHECK ((*itv).first == element1, internal_logic ()); // broken map 1556 const_subiterator_type it ((*itv).second.find (element2)); 1557 if (it == (*itv).second.end ()) 1558 return zero_; 1559 BOOST_UBLAS_CHECK ((*itv).first == element1, internal_logic ()); // broken map 1560 return (*it).second; 1561 } 1562 BOOST_UBLAS_INLINE operator ()(size_type i,size_type j)1563 reference operator () (size_type i, size_type j) { 1564 #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE 1565 const size_type element1 = layout_type::index_M (i, j); 1566 const size_type element2 = layout_type::index_m (i, j); 1567 vector_data_value_type& vd (data () [element1]); 1568 std::pair<subiterator_type, bool> ii (vd.insert (typename array_type::value_type::second_type::value_type (element2, value_type/*zero*/()))); 1569 BOOST_UBLAS_CHECK ((ii.first)->first == element2, internal_logic ()); // broken map 1570 return (ii.first)->second; 1571 #else 1572 return reference (*this, i, j); 1573 #endif 1574 } 1575 1576 // Element assignment 1577 BOOST_UBLAS_INLINE insert_element(size_type i,size_type j,const_reference t)1578 true_reference insert_element (size_type i, size_type j, const_reference t) { 1579 BOOST_UBLAS_CHECK (!find_element (i, j), bad_index ()); // duplicate element 1580 const size_type element1 = layout_type::index_M (i, j); 1581 const size_type element2 = layout_type::index_m (i, j); 1582 1583 vector_data_value_type& vd (data () [element1]); 1584 std::pair<subiterator_type, bool> ii (vd.insert (typename vector_data_value_type::value_type (element2, t))); 1585 BOOST_UBLAS_CHECK ((ii.first)->first == element2, internal_logic ()); // broken map 1586 if (!ii.second) // existing element 1587 (ii.first)->second = t; 1588 return (ii.first)->second; 1589 } 1590 BOOST_UBLAS_INLINE erase_element(size_type i,size_type j)1591 void erase_element (size_type i, size_type j) { 1592 vector_subiterator_type itv (data ().find (layout_type::index_M (i, j))); 1593 if (itv == data ().end ()) 1594 return; 1595 subiterator_type it ((*itv).second.find (layout_type::index_m (i, j))); 1596 if (it == (*itv).second.end ()) 1597 return; 1598 (*itv).second.erase (it); 1599 } 1600 1601 // Zeroing 1602 BOOST_UBLAS_INLINE clear()1603 void clear () { 1604 data ().clear (); 1605 data_ [layout_type::size_M (size1_, size2_)] = vector_data_value_type (); 1606 } 1607 1608 // Assignment 1609 BOOST_UBLAS_INLINE operator =(const mapped_vector_of_mapped_vector & m)1610 mapped_vector_of_mapped_vector &operator = (const mapped_vector_of_mapped_vector &m) { 1611 if (this != &m) { 1612 size1_ = m.size1_; 1613 size2_ = m.size2_; 1614 data () = m.data (); 1615 } 1616 return *this; 1617 } 1618 template<class C> // Container assignment without temporary 1619 BOOST_UBLAS_INLINE operator =(const matrix_container<C> & m)1620 mapped_vector_of_mapped_vector &operator = (const matrix_container<C> &m) { 1621 resize (m ().size1 (), m ().size2 (), false); 1622 assign (m); 1623 return *this; 1624 } 1625 BOOST_UBLAS_INLINE assign_temporary(mapped_vector_of_mapped_vector & m)1626 mapped_vector_of_mapped_vector &assign_temporary (mapped_vector_of_mapped_vector &m) { 1627 swap (m); 1628 return *this; 1629 } 1630 template<class AE> 1631 BOOST_UBLAS_INLINE operator =(const matrix_expression<AE> & ae)1632 mapped_vector_of_mapped_vector &operator = (const matrix_expression<AE> &ae) { 1633 self_type temporary (ae); 1634 return assign_temporary (temporary); 1635 } 1636 template<class AE> 1637 BOOST_UBLAS_INLINE assign(const matrix_expression<AE> & ae)1638 mapped_vector_of_mapped_vector &assign (const matrix_expression<AE> &ae) { 1639 matrix_assign<scalar_assign> (*this, ae); 1640 return *this; 1641 } 1642 template<class AE> 1643 BOOST_UBLAS_INLINE operator +=(const matrix_expression<AE> & ae)1644 mapped_vector_of_mapped_vector& operator += (const matrix_expression<AE> &ae) { 1645 self_type temporary (*this + ae); 1646 return assign_temporary (temporary); 1647 } 1648 template<class C> // Container assignment without temporary 1649 BOOST_UBLAS_INLINE operator +=(const matrix_container<C> & m)1650 mapped_vector_of_mapped_vector &operator += (const matrix_container<C> &m) { 1651 plus_assign (m); 1652 return *this; 1653 } 1654 template<class AE> 1655 BOOST_UBLAS_INLINE plus_assign(const matrix_expression<AE> & ae)1656 mapped_vector_of_mapped_vector &plus_assign (const matrix_expression<AE> &ae) { 1657 matrix_assign<scalar_plus_assign> (*this, ae); 1658 return *this; 1659 } 1660 template<class AE> 1661 BOOST_UBLAS_INLINE operator -=(const matrix_expression<AE> & ae)1662 mapped_vector_of_mapped_vector& operator -= (const matrix_expression<AE> &ae) { 1663 self_type temporary (*this - ae); 1664 return assign_temporary (temporary); 1665 } 1666 template<class C> // Container assignment without temporary 1667 BOOST_UBLAS_INLINE operator -=(const matrix_container<C> & m)1668 mapped_vector_of_mapped_vector &operator -= (const matrix_container<C> &m) { 1669 minus_assign (m); 1670 return *this; 1671 } 1672 template<class AE> 1673 BOOST_UBLAS_INLINE minus_assign(const matrix_expression<AE> & ae)1674 mapped_vector_of_mapped_vector &minus_assign (const matrix_expression<AE> &ae) { 1675 matrix_assign<scalar_minus_assign> (*this, ae); 1676 return *this; 1677 } 1678 template<class AT> 1679 BOOST_UBLAS_INLINE operator *=(const AT & at)1680 mapped_vector_of_mapped_vector& operator *= (const AT &at) { 1681 matrix_assign_scalar<scalar_multiplies_assign> (*this, at); 1682 return *this; 1683 } 1684 template<class AT> 1685 BOOST_UBLAS_INLINE operator /=(const AT & at)1686 mapped_vector_of_mapped_vector& operator /= (const AT &at) { 1687 matrix_assign_scalar<scalar_divides_assign> (*this, at); 1688 return *this; 1689 } 1690 1691 // Swapping 1692 BOOST_UBLAS_INLINE swap(mapped_vector_of_mapped_vector & m)1693 void swap (mapped_vector_of_mapped_vector &m) { 1694 if (this != &m) { 1695 std::swap (size1_, m.size1_); 1696 std::swap (size2_, m.size2_); 1697 data ().swap (m.data ()); 1698 } 1699 } 1700 BOOST_UBLAS_INLINE swap(mapped_vector_of_mapped_vector & m1,mapped_vector_of_mapped_vector & m2)1701 friend void swap (mapped_vector_of_mapped_vector &m1, mapped_vector_of_mapped_vector &m2) { 1702 m1.swap (m2); 1703 } 1704 1705 // Iterator types 1706 private: 1707 // Use storage iterators 1708 typedef typename A::const_iterator vector_const_subiterator_type; 1709 typedef typename A::iterator vector_subiterator_type; 1710 typedef typename A::value_type::second_type::const_iterator const_subiterator_type; 1711 typedef typename A::value_type::second_type::iterator subiterator_type; 1712 1713 BOOST_UBLAS_INLINE at_element(size_type i,size_type j)1714 true_reference at_element (size_type i, size_type j) { 1715 const size_type element1 = layout_type::index_M (i, j); 1716 const size_type element2 = layout_type::index_m (i, j); 1717 vector_subiterator_type itv (data ().find (element1)); 1718 BOOST_UBLAS_CHECK (itv != data ().end(), bad_index ()); 1719 BOOST_UBLAS_CHECK ((*itv).first == element1, internal_logic ()); // broken map 1720 subiterator_type it ((*itv).second.find (element2)); 1721 BOOST_UBLAS_CHECK (it != (*itv).second.end (), bad_index ()); 1722 BOOST_UBLAS_CHECK ((*it).first == element2, internal_logic ()); // broken map 1723 1724 return it->second; 1725 } 1726 1727 public: 1728 class const_iterator1; 1729 class iterator1; 1730 class const_iterator2; 1731 class iterator2; 1732 typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1; 1733 typedef reverse_iterator_base1<iterator1> reverse_iterator1; 1734 typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2; 1735 typedef reverse_iterator_base2<iterator2> reverse_iterator2; 1736 1737 // Element lookup 1738 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. find1(int rank,size_type i,size_type j,int direction=1) const1739 const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const { 1740 BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ()); 1741 for (;;) { 1742 vector_const_subiterator_type itv (data ().lower_bound (layout_type::index_M (i, j))); 1743 vector_const_subiterator_type itv_end (data ().end ()); 1744 if (itv == itv_end) 1745 return const_iterator1 (*this, rank, i, j, itv_end, (*(-- itv)).second.end ()); 1746 1747 const_subiterator_type it ((*itv).second.lower_bound (layout_type::index_m (i, j))); 1748 const_subiterator_type it_end ((*itv).second.end ()); 1749 if (rank == 0) { 1750 // advance to the first available major index 1751 size_type M = itv->first; 1752 size_type m; 1753 if (it != it_end) { 1754 m = it->first; 1755 } else { 1756 m = layout_type::size_m(size1_, size2_); 1757 } 1758 size_type first_i = layout_type::index_M(M,m); 1759 return const_iterator1 (*this, rank, first_i, j, itv, it); 1760 } 1761 if (it != it_end && (*it).first == layout_type::index_m (i, j)) 1762 return const_iterator1 (*this, rank, i, j, itv, it); 1763 if (direction > 0) { 1764 if (layout_type::fast_i ()) { 1765 if (it == it_end) 1766 return const_iterator1 (*this, rank, i, j, itv, it); 1767 i = (*it).first; 1768 } else { 1769 if (i >= size1_) 1770 return const_iterator1 (*this, rank, i, j, itv, it); 1771 ++ i; 1772 } 1773 } else /* if (direction < 0) */ { 1774 if (layout_type::fast_i ()) { 1775 if (it == (*itv).second.begin ()) 1776 return const_iterator1 (*this, rank, i, j, itv, it); 1777 -- it; 1778 i = (*it).first; 1779 } else { 1780 if (i == 0) 1781 return const_iterator1 (*this, rank, i, j, itv, it); 1782 -- i; 1783 } 1784 } 1785 } 1786 } 1787 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. find1(int rank,size_type i,size_type j,int direction=1)1788 iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) { 1789 BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ()); 1790 for (;;) { 1791 vector_subiterator_type itv (data ().lower_bound (layout_type::index_M (i, j))); 1792 vector_subiterator_type itv_end (data ().end ()); 1793 if (itv == itv_end) 1794 return iterator1 (*this, rank, i, j, itv_end, (*(-- itv)).second.end ()); 1795 1796 subiterator_type it ((*itv).second.lower_bound (layout_type::index_m (i, j))); 1797 subiterator_type it_end ((*itv).second.end ()); 1798 if (rank == 0) { 1799 // advance to the first available major index 1800 size_type M = itv->first; 1801 size_type m; 1802 if (it != it_end) { 1803 m = it->first; 1804 } else { 1805 m = layout_type::size_m(size1_, size2_); 1806 } 1807 size_type first_i = layout_type::index_M(M,m); 1808 return iterator1 (*this, rank, first_i, j, itv, it); 1809 } 1810 if (it != it_end && (*it).first == layout_type::index_m (i, j)) 1811 return iterator1 (*this, rank, i, j, itv, it); 1812 if (direction > 0) { 1813 if (layout_type::fast_i ()) { 1814 if (it == it_end) 1815 return iterator1 (*this, rank, i, j, itv, it); 1816 i = (*it).first; 1817 } else { 1818 if (i >= size1_) 1819 return iterator1 (*this, rank, i, j, itv, it); 1820 ++ i; 1821 } 1822 } else /* if (direction < 0) */ { 1823 if (layout_type::fast_i ()) { 1824 if (it == (*itv).second.begin ()) 1825 return iterator1 (*this, rank, i, j, itv, it); 1826 -- it; 1827 i = (*it).first; 1828 } else { 1829 if (i == 0) 1830 return iterator1 (*this, rank, i, j, itv, it); 1831 -- i; 1832 } 1833 } 1834 } 1835 } 1836 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. find2(int rank,size_type i,size_type j,int direction=1) const1837 const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const { 1838 BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ()); 1839 for (;;) { 1840 vector_const_subiterator_type itv (data ().lower_bound (layout_type::index_M (i, j))); 1841 vector_const_subiterator_type itv_end (data ().end ()); 1842 if (itv == itv_end) 1843 return const_iterator2 (*this, rank, i, j, itv_end, (*(-- itv)).second.end ()); 1844 1845 const_subiterator_type it ((*itv).second.lower_bound (layout_type::index_m (i, j))); 1846 const_subiterator_type it_end ((*itv).second.end ()); 1847 if (rank == 0) { 1848 // advance to the first available major index 1849 size_type M = itv->first; 1850 size_type m; 1851 if (it != it_end) { 1852 m = it->first; 1853 } else { 1854 m = layout_type::size_m(size1_, size2_); 1855 } 1856 size_type first_j = layout_type::index_m(M,m); 1857 return const_iterator2 (*this, rank, i, first_j, itv, it); 1858 } 1859 if (it != it_end && (*it).first == layout_type::index_m (i, j)) 1860 return const_iterator2 (*this, rank, i, j, itv, it); 1861 if (direction > 0) { 1862 if (layout_type::fast_j ()) { 1863 if (it == it_end) 1864 return const_iterator2 (*this, rank, i, j, itv, it); 1865 j = (*it).first; 1866 } else { 1867 if (j >= size2_) 1868 return const_iterator2 (*this, rank, i, j, itv, it); 1869 ++ j; 1870 } 1871 } else /* if (direction < 0) */ { 1872 if (layout_type::fast_j ()) { 1873 if (it == (*itv).second.begin ()) 1874 return const_iterator2 (*this, rank, i, j, itv, it); 1875 -- it; 1876 j = (*it).first; 1877 } else { 1878 if (j == 0) 1879 return const_iterator2 (*this, rank, i, j, itv, it); 1880 -- j; 1881 } 1882 } 1883 } 1884 } 1885 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. find2(int rank,size_type i,size_type j,int direction=1)1886 iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) { 1887 BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ()); 1888 for (;;) { 1889 vector_subiterator_type itv (data ().lower_bound (layout_type::index_M (i, j))); 1890 vector_subiterator_type itv_end (data ().end ()); 1891 if (itv == itv_end) 1892 return iterator2 (*this, rank, i, j, itv_end, (*(-- itv)).second.end ()); 1893 1894 subiterator_type it ((*itv).second.lower_bound (layout_type::index_m (i, j))); 1895 subiterator_type it_end ((*itv).second.end ()); 1896 if (rank == 0) { 1897 // advance to the first available major index 1898 size_type M = itv->first; 1899 size_type m; 1900 if (it != it_end) { 1901 m = it->first; 1902 } else { 1903 m = layout_type::size_m(size1_, size2_); 1904 } 1905 size_type first_j = layout_type::index_m(M,m); 1906 return iterator2 (*this, rank, i, first_j, itv, it); 1907 } 1908 if (it != it_end && (*it).first == layout_type::index_m (i, j)) 1909 return iterator2 (*this, rank, i, j, itv, it); 1910 if (direction > 0) { 1911 if (layout_type::fast_j ()) { 1912 if (it == it_end) 1913 return iterator2 (*this, rank, i, j, itv, it); 1914 j = (*it).first; 1915 } else { 1916 if (j >= size2_) 1917 return iterator2 (*this, rank, i, j, itv, it); 1918 ++ j; 1919 } 1920 } else /* if (direction < 0) */ { 1921 if (layout_type::fast_j ()) { 1922 if (it == (*itv).second.begin ()) 1923 return iterator2 (*this, rank, i, j, itv, it); 1924 -- it; 1925 j = (*it).first; 1926 } else { 1927 if (j == 0) 1928 return iterator2 (*this, rank, i, j, itv, it); 1929 -- j; 1930 } 1931 } 1932 } 1933 } 1934 1935 class const_iterator1: 1936 public container_const_reference<mapped_vector_of_mapped_vector>, 1937 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag, 1938 const_iterator1, value_type> { 1939 public: 1940 typedef typename mapped_vector_of_mapped_vector::value_type value_type; 1941 typedef typename mapped_vector_of_mapped_vector::difference_type difference_type; 1942 typedef typename mapped_vector_of_mapped_vector::const_reference reference; 1943 typedef const typename mapped_vector_of_mapped_vector::pointer pointer; 1944 1945 typedef const_iterator2 dual_iterator_type; 1946 typedef const_reverse_iterator2 dual_reverse_iterator_type; 1947 1948 // Construction and destruction 1949 BOOST_UBLAS_INLINE const_iterator1()1950 const_iterator1 (): 1951 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {} 1952 BOOST_UBLAS_INLINE const_iterator1(const self_type & m,int rank,size_type i,size_type j,const vector_const_subiterator_type & itv,const const_subiterator_type & it)1953 const_iterator1 (const self_type &m, int rank, size_type i, size_type j, const vector_const_subiterator_type &itv, const const_subiterator_type &it): 1954 container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {} 1955 BOOST_UBLAS_INLINE const_iterator1(const iterator1 & it)1956 const_iterator1 (const iterator1 &it): 1957 container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {} 1958 1959 // Arithmetic 1960 BOOST_UBLAS_INLINE operator ++()1961 const_iterator1 &operator ++ () { 1962 if (rank_ == 1 && layout_type::fast_i ()) 1963 ++ it_; 1964 else { 1965 const self_type &m = (*this) (); 1966 if (rank_ == 0) { 1967 ++ itv_; 1968 i_ = itv_->first; 1969 } else { 1970 i_ = index1 () + 1; 1971 } 1972 if (rank_ == 1 && ++ itv_ == m.end1 ().itv_) 1973 *this = m.find1 (rank_, i_, j_, 1); 1974 else if (rank_ == 1) { 1975 it_ = (*itv_).second.begin (); 1976 if (it_ == (*itv_).second.end () || index2 () != j_) 1977 *this = m.find1 (rank_, i_, j_, 1); 1978 } 1979 } 1980 return *this; 1981 } 1982 BOOST_UBLAS_INLINE operator --()1983 const_iterator1 &operator -- () { 1984 if (rank_ == 1 && layout_type::fast_i ()) 1985 -- it_; 1986 else { 1987 const self_type &m = (*this) (); 1988 if (rank_ == 0) { 1989 -- itv_; 1990 i_ = itv_->first; 1991 } else { 1992 i_ = index1 () - 1; 1993 } 1994 // FIXME: this expression should never become true! 1995 if (rank_ == 1 && -- itv_ == m.end1 ().itv_) 1996 *this = m.find1 (rank_, i_, j_, -1); 1997 else if (rank_ == 1) { 1998 it_ = (*itv_).second.begin (); 1999 if (it_ == (*itv_).second.end () || index2 () != j_) 2000 *this = m.find1 (rank_, i_, j_, -1); 2001 } 2002 } 2003 return *this; 2004 } 2005 2006 // Dereference 2007 BOOST_UBLAS_INLINE operator *() const2008 const_reference operator * () const { 2009 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ()); 2010 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ()); 2011 if (rank_ == 1) { 2012 return (*it_).second; 2013 } else { 2014 return (*this) () (i_, j_); 2015 } 2016 } 2017 2018 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 2019 BOOST_UBLAS_INLINE 2020 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2021 typename self_type:: 2022 #endif begin() const2023 const_iterator2 begin () const { 2024 const self_type &m = (*this) (); 2025 return m.find2 (1, index1 (), 0); 2026 } 2027 BOOST_UBLAS_INLINE 2028 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2029 typename self_type:: 2030 #endif cbegin() const2031 const_iterator2 cbegin () const { 2032 return begin (); 2033 } 2034 BOOST_UBLAS_INLINE 2035 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2036 typename self_type:: 2037 #endif end() const2038 const_iterator2 end () const { 2039 const self_type &m = (*this) (); 2040 return m.find2 (1, index1 (), m.size2 ()); 2041 } 2042 BOOST_UBLAS_INLINE 2043 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2044 typename self_type:: 2045 #endif cend() const2046 const_iterator2 cend () const { 2047 return end (); 2048 } 2049 BOOST_UBLAS_INLINE 2050 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2051 typename self_type:: 2052 #endif rbegin() const2053 const_reverse_iterator2 rbegin () const { 2054 return const_reverse_iterator2 (end ()); 2055 } 2056 BOOST_UBLAS_INLINE 2057 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2058 typename self_type:: 2059 #endif crbegin() const2060 const_reverse_iterator2 crbegin () const { 2061 return rbegin (); 2062 } 2063 BOOST_UBLAS_INLINE 2064 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2065 typename self_type:: 2066 #endif rend() const2067 const_reverse_iterator2 rend () const { 2068 return const_reverse_iterator2 (begin ()); 2069 } 2070 BOOST_UBLAS_INLINE 2071 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2072 typename self_type:: 2073 #endif crend() const2074 const_reverse_iterator2 crend () const { 2075 return rend (); 2076 } 2077 #endif 2078 2079 // Indices 2080 BOOST_UBLAS_INLINE index1() const2081 size_type index1 () const { 2082 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ()); 2083 if (rank_ == 1) { 2084 BOOST_UBLAS_CHECK (layout_type::index_M ((*itv_).first, (*it_).first) < (*this) ().size1 (), bad_index ()); 2085 return layout_type::index_M ((*itv_).first, (*it_).first); 2086 } else { 2087 return i_; 2088 } 2089 } 2090 BOOST_UBLAS_INLINE index2() const2091 size_type index2 () const { 2092 if (rank_ == 1) { 2093 BOOST_UBLAS_CHECK (layout_type::index_m ((*itv_).first, (*it_).first) < (*this) ().size2 (), bad_index ()); 2094 return layout_type::index_m ((*itv_).first, (*it_).first); 2095 } else { 2096 return j_; 2097 } 2098 } 2099 2100 // Assignment 2101 BOOST_UBLAS_INLINE operator =(const const_iterator1 & it)2102 const_iterator1 &operator = (const const_iterator1 &it) { 2103 container_const_reference<self_type>::assign (&it ()); 2104 rank_ = it.rank_; 2105 i_ = it.i_; 2106 j_ = it.j_; 2107 itv_ = it.itv_; 2108 it_ = it.it_; 2109 return *this; 2110 } 2111 2112 // Comparison 2113 BOOST_UBLAS_INLINE operator ==(const const_iterator1 & it) const2114 bool operator == (const const_iterator1 &it) const { 2115 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 2116 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ()); 2117 if (rank_ == 1 || it.rank_ == 1) { 2118 return it_ == it.it_; 2119 } else { 2120 return i_ == it.i_ && j_ == it.j_; 2121 } 2122 } 2123 2124 private: 2125 int rank_; 2126 size_type i_; 2127 size_type j_; 2128 vector_const_subiterator_type itv_; 2129 const_subiterator_type it_; 2130 }; 2131 2132 BOOST_UBLAS_INLINE begin1() const2133 const_iterator1 begin1 () const { 2134 return find1 (0, 0, 0); 2135 } 2136 BOOST_UBLAS_INLINE cbegin1() const2137 const_iterator1 cbegin1 () const { 2138 return begin1 (); 2139 } 2140 BOOST_UBLAS_INLINE end1() const2141 const_iterator1 end1 () const { 2142 return find1 (0, size1_, 0); 2143 } 2144 BOOST_UBLAS_INLINE cend1() const2145 const_iterator1 cend1 () const { 2146 return end1 (); 2147 } 2148 2149 class iterator1: 2150 public container_reference<mapped_vector_of_mapped_vector>, 2151 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag, 2152 iterator1, value_type> { 2153 public: 2154 typedef typename mapped_vector_of_mapped_vector::value_type value_type; 2155 typedef typename mapped_vector_of_mapped_vector::difference_type difference_type; 2156 typedef typename mapped_vector_of_mapped_vector::true_reference reference; 2157 typedef typename mapped_vector_of_mapped_vector::pointer pointer; 2158 2159 typedef iterator2 dual_iterator_type; 2160 typedef reverse_iterator2 dual_reverse_iterator_type; 2161 2162 // Construction and destruction 2163 BOOST_UBLAS_INLINE iterator1()2164 iterator1 (): 2165 container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {} 2166 BOOST_UBLAS_INLINE iterator1(self_type & m,int rank,size_type i,size_type j,const vector_subiterator_type & itv,const subiterator_type & it)2167 iterator1 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it): 2168 container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {} 2169 2170 // Arithmetic 2171 BOOST_UBLAS_INLINE operator ++()2172 iterator1 &operator ++ () { 2173 if (rank_ == 1 && layout_type::fast_i ()) 2174 ++ it_; 2175 else { 2176 self_type &m = (*this) (); 2177 if (rank_ == 0) { 2178 ++ itv_; 2179 i_ = itv_->first; 2180 } else { 2181 i_ = index1 () + 1; 2182 } 2183 if (rank_ == 1 && ++ itv_ == m.end1 ().itv_) 2184 *this = m.find1 (rank_, i_, j_, 1); 2185 else if (rank_ == 1) { 2186 it_ = (*itv_).second.begin (); 2187 if (it_ == (*itv_).second.end () || index2 () != j_) 2188 *this = m.find1 (rank_, i_, j_, 1); 2189 } 2190 } 2191 return *this; 2192 } 2193 BOOST_UBLAS_INLINE operator --()2194 iterator1 &operator -- () { 2195 if (rank_ == 1 && layout_type::fast_i ()) 2196 -- it_; 2197 else { 2198 self_type &m = (*this) (); 2199 if (rank_ == 0) { 2200 -- itv_; 2201 i_ = itv_->first; 2202 } else { 2203 i_ = index1 () - 1; 2204 } 2205 // FIXME: this expression should never become true! 2206 if (rank_ == 1 && -- itv_ == m.end1 ().itv_) 2207 *this = m.find1 (rank_, i_, j_, -1); 2208 else if (rank_ == 1) { 2209 it_ = (*itv_).second.begin (); 2210 if (it_ == (*itv_).second.end () || index2 () != j_) 2211 *this = m.find1 (rank_, i_, j_, -1); 2212 } 2213 } 2214 return *this; 2215 } 2216 2217 // Dereference 2218 BOOST_UBLAS_INLINE operator *() const2219 reference operator * () const { 2220 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ()); 2221 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ()); 2222 if (rank_ == 1) { 2223 return (*it_).second; 2224 } else { 2225 return (*this) ().at_element (i_, j_); 2226 } 2227 } 2228 2229 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 2230 BOOST_UBLAS_INLINE 2231 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2232 typename self_type:: 2233 #endif begin() const2234 iterator2 begin () const { 2235 self_type &m = (*this) (); 2236 return m.find2 (1, index1 (), 0); 2237 } 2238 BOOST_UBLAS_INLINE 2239 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2240 typename self_type:: 2241 #endif end() const2242 iterator2 end () const { 2243 self_type &m = (*this) (); 2244 return m.find2 (1, index1 (), m.size2 ()); 2245 } 2246 BOOST_UBLAS_INLINE 2247 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2248 typename self_type:: 2249 #endif rbegin() const2250 reverse_iterator2 rbegin () const { 2251 return reverse_iterator2 (end ()); 2252 } 2253 BOOST_UBLAS_INLINE 2254 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2255 typename self_type:: 2256 #endif rend() const2257 reverse_iterator2 rend () const { 2258 return reverse_iterator2 (begin ()); 2259 } 2260 #endif 2261 2262 // Indices 2263 BOOST_UBLAS_INLINE index1() const2264 size_type index1 () const { 2265 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ()); 2266 if (rank_ == 1) { 2267 BOOST_UBLAS_CHECK (layout_type::index_M ((*itv_).first, (*it_).first) < (*this) ().size1 (), bad_index ()); 2268 return layout_type::index_M ((*itv_).first, (*it_).first); 2269 } else { 2270 return i_; 2271 } 2272 } 2273 BOOST_UBLAS_INLINE index2() const2274 size_type index2 () const { 2275 if (rank_ == 1) { 2276 BOOST_UBLAS_CHECK (layout_type::index_m ((*itv_).first, (*it_).first) < (*this) ().size2 (), bad_index ()); 2277 return layout_type::index_m ((*itv_).first, (*it_).first); 2278 } else { 2279 return j_; 2280 } 2281 } 2282 2283 // Assignment 2284 BOOST_UBLAS_INLINE operator =(const iterator1 & it)2285 iterator1 &operator = (const iterator1 &it) { 2286 container_reference<self_type>::assign (&it ()); 2287 rank_ = it.rank_; 2288 i_ = it.i_; 2289 j_ = it.j_; 2290 itv_ = it.itv_; 2291 it_ = it.it_; 2292 return *this; 2293 } 2294 2295 // Comparison 2296 BOOST_UBLAS_INLINE operator ==(const iterator1 & it) const2297 bool operator == (const iterator1 &it) const { 2298 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 2299 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ()); 2300 if (rank_ == 1 || it.rank_ == 1) { 2301 return it_ == it.it_; 2302 } else { 2303 return i_ == it.i_ && j_ == it.j_; 2304 } 2305 } 2306 2307 private: 2308 int rank_; 2309 size_type i_; 2310 size_type j_; 2311 vector_subiterator_type itv_; 2312 subiterator_type it_; 2313 2314 friend class const_iterator1; 2315 }; 2316 2317 BOOST_UBLAS_INLINE begin1()2318 iterator1 begin1 () { 2319 return find1 (0, 0, 0); 2320 } 2321 BOOST_UBLAS_INLINE end1()2322 iterator1 end1 () { 2323 return find1 (0, size1_, 0); 2324 } 2325 2326 class const_iterator2: 2327 public container_const_reference<mapped_vector_of_mapped_vector>, 2328 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag, 2329 const_iterator2, value_type> { 2330 public: 2331 typedef typename mapped_vector_of_mapped_vector::value_type value_type; 2332 typedef typename mapped_vector_of_mapped_vector::difference_type difference_type; 2333 typedef typename mapped_vector_of_mapped_vector::const_reference reference; 2334 typedef const typename mapped_vector_of_mapped_vector::pointer pointer; 2335 2336 typedef const_iterator1 dual_iterator_type; 2337 typedef const_reverse_iterator1 dual_reverse_iterator_type; 2338 2339 // Construction and destruction 2340 BOOST_UBLAS_INLINE const_iterator2()2341 const_iterator2 (): 2342 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {} 2343 BOOST_UBLAS_INLINE const_iterator2(const self_type & m,int rank,size_type i,size_type j,const vector_const_subiterator_type & itv,const const_subiterator_type & it)2344 const_iterator2 (const self_type &m, int rank, size_type i, size_type j, const vector_const_subiterator_type &itv, const const_subiterator_type &it): 2345 container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {} 2346 BOOST_UBLAS_INLINE const_iterator2(const iterator2 & it)2347 const_iterator2 (const iterator2 &it): 2348 container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {} 2349 2350 // Arithmetic 2351 BOOST_UBLAS_INLINE operator ++()2352 const_iterator2 &operator ++ () { 2353 if (rank_ == 1 && layout_type::fast_j ()) 2354 ++ it_; 2355 else { 2356 const self_type &m = (*this) (); 2357 if (rank_ == 0) { 2358 ++ itv_; 2359 j_ = itv_->first; 2360 } else { 2361 j_ = index2 () + 1; 2362 } 2363 if (rank_ == 1 && ++ itv_ == m.end2 ().itv_) 2364 *this = m.find2 (rank_, i_, j_, 1); 2365 else if (rank_ == 1) { 2366 it_ = (*itv_).second.begin (); 2367 if (it_ == (*itv_).second.end () || index1 () != i_) 2368 *this = m.find2 (rank_, i_, j_, 1); 2369 } 2370 } 2371 return *this; 2372 } 2373 BOOST_UBLAS_INLINE operator --()2374 const_iterator2 &operator -- () { 2375 if (rank_ == 1 && layout_type::fast_j ()) 2376 -- it_; 2377 else { 2378 const self_type &m = (*this) (); 2379 if (rank_ == 0) { 2380 -- itv_; 2381 j_ = itv_->first; 2382 } else { 2383 j_ = index2 () - 1; 2384 } 2385 // FIXME: this expression should never become true! 2386 if (rank_ == 1 && -- itv_ == m.end2 ().itv_) 2387 *this = m.find2 (rank_, i_, j_, -1); 2388 else if (rank_ == 1) { 2389 it_ = (*itv_).second.begin (); 2390 if (it_ == (*itv_).second.end () || index1 () != i_) 2391 *this = m.find2 (rank_, i_, j_, -1); 2392 } 2393 } 2394 return *this; 2395 } 2396 2397 // Dereference 2398 BOOST_UBLAS_INLINE operator *() const2399 const_reference operator * () const { 2400 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ()); 2401 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ()); 2402 if (rank_ == 1) { 2403 return (*it_).second; 2404 } else { 2405 return (*this) () (i_, j_); 2406 } 2407 } 2408 2409 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 2410 BOOST_UBLAS_INLINE 2411 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2412 typename self_type:: 2413 #endif begin() const2414 const_iterator1 begin () const { 2415 const self_type &m = (*this) (); 2416 return m.find1 (1, 0, index2 ()); 2417 } 2418 BOOST_UBLAS_INLINE 2419 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2420 typename self_type:: 2421 #endif cbegin() const2422 const_iterator1 cbegin () const { 2423 return begin (); 2424 } 2425 BOOST_UBLAS_INLINE 2426 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2427 typename self_type:: 2428 #endif end() const2429 const_iterator1 end () const { 2430 const self_type &m = (*this) (); 2431 return m.find1 (1, m.size1 (), index2 ()); 2432 } 2433 BOOST_UBLAS_INLINE 2434 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2435 typename self_type:: 2436 #endif cend() const2437 const_iterator1 cend () const { 2438 return end (); 2439 } 2440 BOOST_UBLAS_INLINE 2441 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2442 typename self_type:: 2443 #endif rbegin() const2444 const_reverse_iterator1 rbegin () const { 2445 return const_reverse_iterator1 (end ()); 2446 } 2447 BOOST_UBLAS_INLINE 2448 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2449 typename self_type:: 2450 #endif crbegin() const2451 const_reverse_iterator1 crbegin () const { 2452 return rbegin (); 2453 } 2454 BOOST_UBLAS_INLINE 2455 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2456 typename self_type:: 2457 #endif rend() const2458 const_reverse_iterator1 rend () const { 2459 return const_reverse_iterator1 (begin ()); 2460 } 2461 BOOST_UBLAS_INLINE 2462 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2463 typename self_type:: 2464 #endif crend() const2465 const_reverse_iterator1 crend () const { 2466 return rend (); 2467 } 2468 #endif 2469 2470 // Indices 2471 BOOST_UBLAS_INLINE index1() const2472 size_type index1 () const { 2473 if (rank_ == 1) { 2474 BOOST_UBLAS_CHECK (layout_type::index_M ((*itv_).first, (*it_).first) < (*this) ().size1 (), bad_index ()); 2475 return layout_type::index_M ((*itv_).first, (*it_).first); 2476 } else { 2477 return i_; 2478 } 2479 } 2480 BOOST_UBLAS_INLINE index2() const2481 size_type index2 () const { 2482 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ()); 2483 if (rank_ == 1) { 2484 BOOST_UBLAS_CHECK (layout_type::index_m ((*itv_).first, (*it_).first) < (*this) ().size2 (), bad_index ()); 2485 return layout_type::index_m ((*itv_).first, (*it_).first); 2486 } else { 2487 return j_; 2488 } 2489 } 2490 2491 // Assignment 2492 BOOST_UBLAS_INLINE operator =(const const_iterator2 & it)2493 const_iterator2 &operator = (const const_iterator2 &it) { 2494 container_const_reference<self_type>::assign (&it ()); 2495 rank_ = it.rank_; 2496 i_ = it.i_; 2497 j_ = it.j_; 2498 itv_ = it.itv_; 2499 it_ = it.it_; 2500 return *this; 2501 } 2502 2503 // Comparison 2504 BOOST_UBLAS_INLINE operator ==(const const_iterator2 & it) const2505 bool operator == (const const_iterator2 &it) const { 2506 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 2507 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ()); 2508 if (rank_ == 1 || it.rank_ == 1) { 2509 return it_ == it.it_; 2510 } else { 2511 return i_ == it.i_ && j_ == it.j_; 2512 } 2513 } 2514 2515 private: 2516 int rank_; 2517 size_type i_; 2518 size_type j_; 2519 vector_const_subiterator_type itv_; 2520 const_subiterator_type it_; 2521 }; 2522 2523 BOOST_UBLAS_INLINE begin2() const2524 const_iterator2 begin2 () const { 2525 return find2 (0, 0, 0); 2526 } 2527 BOOST_UBLAS_INLINE cbegin2() const2528 const_iterator2 cbegin2 () const { 2529 return begin2 (); 2530 } 2531 BOOST_UBLAS_INLINE end2() const2532 const_iterator2 end2 () const { 2533 return find2 (0, 0, size2_); 2534 } 2535 BOOST_UBLAS_INLINE cend2() const2536 const_iterator2 cend2 () const { 2537 return end2 (); 2538 } 2539 2540 class iterator2: 2541 public container_reference<mapped_vector_of_mapped_vector>, 2542 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag, 2543 iterator2, value_type> { 2544 public: 2545 typedef typename mapped_vector_of_mapped_vector::value_type value_type; 2546 typedef typename mapped_vector_of_mapped_vector::difference_type difference_type; 2547 typedef typename mapped_vector_of_mapped_vector::true_reference reference; 2548 typedef typename mapped_vector_of_mapped_vector::pointer pointer; 2549 2550 typedef iterator1 dual_iterator_type; 2551 typedef reverse_iterator1 dual_reverse_iterator_type; 2552 2553 // Construction and destruction 2554 BOOST_UBLAS_INLINE iterator2()2555 iterator2 (): 2556 container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {} 2557 BOOST_UBLAS_INLINE iterator2(self_type & m,int rank,size_type i,size_type j,const vector_subiterator_type & itv,const subiterator_type & it)2558 iterator2 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it): 2559 container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {} 2560 2561 // Arithmetic 2562 BOOST_UBLAS_INLINE operator ++()2563 iterator2 &operator ++ () { 2564 if (rank_ == 1 && layout_type::fast_j ()) 2565 ++ it_; 2566 else { 2567 self_type &m = (*this) (); 2568 if (rank_ == 0) { 2569 ++ itv_; 2570 j_ = itv_->first; 2571 } else { 2572 j_ = index2 () + 1; 2573 } 2574 if (rank_ == 1 && ++ itv_ == m.end2 ().itv_) 2575 *this = m.find2 (rank_, i_, j_, 1); 2576 else if (rank_ == 1) { 2577 it_ = (*itv_).second.begin (); 2578 if (it_ == (*itv_).second.end () || index1 () != i_) 2579 *this = m.find2 (rank_, i_, j_, 1); 2580 } 2581 } 2582 return *this; 2583 } 2584 BOOST_UBLAS_INLINE operator --()2585 iterator2 &operator -- () { 2586 if (rank_ == 1 && layout_type::fast_j ()) 2587 -- it_; 2588 else { 2589 self_type &m = (*this) (); 2590 if (rank_ == 0) { 2591 -- itv_; 2592 j_ = itv_->first; 2593 } else { 2594 j_ = index2 () - 1; 2595 } 2596 // FIXME: this expression should never become true! 2597 if (rank_ == 1 && -- itv_ == m.end2 ().itv_) 2598 *this = m.find2 (rank_, i_, j_, -1); 2599 else if (rank_ == 1) { 2600 it_ = (*itv_).second.begin (); 2601 if (it_ == (*itv_).second.end () || index1 () != i_) 2602 *this = m.find2 (rank_, i_, j_, -1); 2603 } 2604 } 2605 return *this; 2606 } 2607 2608 // Dereference 2609 BOOST_UBLAS_INLINE operator *() const2610 reference operator * () const { 2611 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ()); 2612 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ()); 2613 if (rank_ == 1) { 2614 return (*it_).second; 2615 } else { 2616 return (*this) ().at_element (i_, j_); 2617 } 2618 } 2619 2620 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 2621 BOOST_UBLAS_INLINE 2622 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2623 typename self_type:: 2624 #endif begin() const2625 iterator1 begin () const { 2626 self_type &m = (*this) (); 2627 return m.find1 (1, 0, index2 ()); 2628 } 2629 BOOST_UBLAS_INLINE 2630 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2631 typename self_type:: 2632 #endif end() const2633 iterator1 end () const { 2634 self_type &m = (*this) (); 2635 return m.find1 (1, m.size1 (), index2 ()); 2636 } 2637 BOOST_UBLAS_INLINE 2638 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2639 typename self_type:: 2640 #endif rbegin() const2641 reverse_iterator1 rbegin () const { 2642 return reverse_iterator1 (end ()); 2643 } 2644 BOOST_UBLAS_INLINE 2645 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 2646 typename self_type:: 2647 #endif rend() const2648 reverse_iterator1 rend () const { 2649 return reverse_iterator1 (begin ()); 2650 } 2651 #endif 2652 2653 // Indices 2654 BOOST_UBLAS_INLINE index1() const2655 size_type index1 () const { 2656 if (rank_ == 1) { 2657 BOOST_UBLAS_CHECK (layout_type::index_M ((*itv_).first, (*it_).first) < (*this) ().size1 (), bad_index ()); 2658 return layout_type::index_M ((*itv_).first, (*it_).first); 2659 } else { 2660 return i_; 2661 } 2662 } 2663 BOOST_UBLAS_INLINE index2() const2664 size_type index2 () const { 2665 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ()); 2666 if (rank_ == 1) { 2667 BOOST_UBLAS_CHECK (layout_type::index_m ((*itv_).first, (*it_).first) < (*this) ().size2 (), bad_index ()); 2668 return layout_type::index_m ((*itv_).first, (*it_).first); 2669 } else { 2670 return j_; 2671 } 2672 } 2673 2674 // Assignment 2675 BOOST_UBLAS_INLINE operator =(const iterator2 & it)2676 iterator2 &operator = (const iterator2 &it) { 2677 container_reference<self_type>::assign (&it ()); 2678 rank_ = it.rank_; 2679 i_ = it.i_; 2680 j_ = it.j_; 2681 itv_ = it.itv_; 2682 it_ = it.it_; 2683 return *this; 2684 } 2685 2686 // Comparison 2687 BOOST_UBLAS_INLINE operator ==(const iterator2 & it) const2688 bool operator == (const iterator2 &it) const { 2689 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 2690 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ()); 2691 if (rank_ == 1 || it.rank_ == 1) { 2692 return it_ == it.it_; 2693 } else { 2694 return i_ == it.i_ && j_ == it.j_; 2695 } 2696 } 2697 2698 private: 2699 int rank_; 2700 size_type i_; 2701 size_type j_; 2702 vector_subiterator_type itv_; 2703 subiterator_type it_; 2704 2705 friend class const_iterator2; 2706 }; 2707 2708 BOOST_UBLAS_INLINE begin2()2709 iterator2 begin2 () { 2710 return find2 (0, 0, 0); 2711 } 2712 BOOST_UBLAS_INLINE end2()2713 iterator2 end2 () { 2714 return find2 (0, 0, size2_); 2715 } 2716 2717 // Reverse iterators 2718 2719 BOOST_UBLAS_INLINE rbegin1() const2720 const_reverse_iterator1 rbegin1 () const { 2721 return const_reverse_iterator1 (end1 ()); 2722 } 2723 BOOST_UBLAS_INLINE crbegin1() const2724 const_reverse_iterator1 crbegin1 () const { 2725 return rbegin1 (); 2726 } 2727 BOOST_UBLAS_INLINE rend1() const2728 const_reverse_iterator1 rend1 () const { 2729 return const_reverse_iterator1 (begin1 ()); 2730 } 2731 BOOST_UBLAS_INLINE crend1() const2732 const_reverse_iterator1 crend1 () const { 2733 return rend1 (); 2734 } 2735 2736 BOOST_UBLAS_INLINE rbegin1()2737 reverse_iterator1 rbegin1 () { 2738 return reverse_iterator1 (end1 ()); 2739 } 2740 BOOST_UBLAS_INLINE rend1()2741 reverse_iterator1 rend1 () { 2742 return reverse_iterator1 (begin1 ()); 2743 } 2744 2745 BOOST_UBLAS_INLINE rbegin2() const2746 const_reverse_iterator2 rbegin2 () const { 2747 return const_reverse_iterator2 (end2 ()); 2748 } 2749 BOOST_UBLAS_INLINE crbegin2() const2750 const_reverse_iterator2 crbegin2 () const { 2751 return rbegin2 (); 2752 } 2753 BOOST_UBLAS_INLINE rend2() const2754 const_reverse_iterator2 rend2 () const { 2755 return const_reverse_iterator2 (begin2 ()); 2756 } 2757 BOOST_UBLAS_INLINE crend2() const2758 const_reverse_iterator2 crend2 () const { 2759 return rend2 (); 2760 } 2761 2762 BOOST_UBLAS_INLINE rbegin2()2763 reverse_iterator2 rbegin2 () { 2764 return reverse_iterator2 (end2 ()); 2765 } 2766 BOOST_UBLAS_INLINE rend2()2767 reverse_iterator2 rend2 () { 2768 return reverse_iterator2 (begin2 ()); 2769 } 2770 2771 // Serialization 2772 template<class Archive> serialize(Archive & ar,const unsigned int)2773 void serialize(Archive & ar, const unsigned int /* file_version */){ 2774 serialization::collection_size_type s1 (size1_); 2775 serialization::collection_size_type s2 (size2_); 2776 ar & serialization::make_nvp("size1",s1); 2777 ar & serialization::make_nvp("size2",s2); 2778 if (Archive::is_loading::value) { 2779 size1_ = s1; 2780 size2_ = s2; 2781 } 2782 ar & serialization::make_nvp("data", data_); 2783 } 2784 2785 private: 2786 size_type size1_; 2787 size_type size2_; 2788 array_type data_; 2789 static const value_type zero_; 2790 }; 2791 2792 template<class T, class L, class A> 2793 const typename mapped_vector_of_mapped_vector<T, L, A>::value_type mapped_vector_of_mapped_vector<T, L, A>::zero_ = value_type/*zero*/(); 2794 2795 2796 // Comperssed array based sparse matrix class 2797 // Thanks to Kresimir Fresl for extending this to cover different index bases. 2798 template<class T, class L, std::size_t IB, class IA, class TA> 2799 class compressed_matrix: 2800 public matrix_container<compressed_matrix<T, L, IB, IA, TA> > { 2801 2802 typedef T &true_reference; 2803 typedef T *pointer; 2804 typedef const T *const_pointer; 2805 typedef L layout_type; 2806 typedef compressed_matrix<T, L, IB, IA, TA> self_type; 2807 public: 2808 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS 2809 using matrix_container<self_type>::operator (); 2810 #endif 2811 // ISSUE require type consistency check 2812 // is_convertable (IA::size_type, TA::size_type) 2813 typedef typename IA::value_type size_type; 2814 // size_type for the data arrays. 2815 typedef typename IA::size_type array_size_type; 2816 // FIXME difference type for sparse storage iterators should it be in the container? 2817 typedef typename IA::difference_type difference_type; 2818 typedef T value_type; 2819 typedef const T &const_reference; 2820 #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE 2821 typedef T &reference; 2822 #else 2823 typedef sparse_matrix_element<self_type> reference; 2824 #endif 2825 typedef IA index_array_type; 2826 typedef TA value_array_type; 2827 typedef const matrix_reference<const self_type> const_closure_type; 2828 typedef matrix_reference<self_type> closure_type; 2829 typedef compressed_vector<T, IB, IA, TA> vector_temporary_type; 2830 typedef self_type matrix_temporary_type; 2831 typedef sparse_tag storage_category; 2832 typedef typename L::orientation_category orientation_category; 2833 2834 // Construction and destruction 2835 BOOST_UBLAS_INLINE compressed_matrix()2836 compressed_matrix (): 2837 matrix_container<self_type> (), 2838 size1_ (0), size2_ (0), capacity_ (restrict_capacity (0)), 2839 filled1_ (1), filled2_ (0), 2840 index1_data_ (layout_type::size_M (size1_, size2_) + 1), index2_data_ (capacity_), value_data_ (capacity_) { 2841 index1_data_ [filled1_ - 1] = k_based (filled2_); 2842 storage_invariants (); 2843 } 2844 BOOST_UBLAS_INLINE compressed_matrix(size_type size1,size_type size2,size_type non_zeros=0)2845 compressed_matrix (size_type size1, size_type size2, size_type non_zeros = 0): 2846 matrix_container<self_type> (), 2847 size1_ (size1), size2_ (size2), capacity_ (restrict_capacity (non_zeros)), 2848 filled1_ (1), filled2_ (0), 2849 index1_data_ (layout_type::size_M (size1_, size2_) + 1), index2_data_ (capacity_), value_data_ (capacity_) { 2850 index1_data_ [filled1_ - 1] = k_based (filled2_); 2851 storage_invariants (); 2852 } 2853 BOOST_UBLAS_INLINE compressed_matrix(const compressed_matrix & m)2854 compressed_matrix (const compressed_matrix &m): 2855 matrix_container<self_type> (), 2856 size1_ (m.size1_), size2_ (m.size2_), capacity_ (m.capacity_), 2857 filled1_ (m.filled1_), filled2_ (m.filled2_), 2858 index1_data_ (m.index1_data_), index2_data_ (m.index2_data_), value_data_ (m.value_data_) { 2859 storage_invariants (); 2860 } 2861 2862 BOOST_UBLAS_INLINE compressed_matrix(const coordinate_matrix<T,L,IB,IA,TA> & m)2863 compressed_matrix (const coordinate_matrix<T, L, IB, IA, TA> &m): 2864 matrix_container<self_type> (), 2865 size1_ (m.size1()), size2_ (m.size2()), 2866 index1_data_ (layout_type::size_M (size1_, size2_) + 1) 2867 { 2868 m.sort(); 2869 reserve(m.nnz(), false); 2870 filled2_ = m.nnz(); 2871 const_subiterator_type i_start = m.index1_data().begin(); 2872 const_subiterator_type i_end = (i_start + filled2_); 2873 const_subiterator_type i = i_start; 2874 size_type r = 1; 2875 for (; (r < layout_type::size_M (size1_, size2_)) && (i != i_end); ++r) { 2876 i = std::lower_bound(i, i_end, r); 2877 index1_data_[r] = k_based( i - i_start ); 2878 } 2879 filled1_ = r + 1; 2880 std::copy( m.index2_data().begin(), m.index2_data().begin() + filled2_, index2_data_.begin()); 2881 std::copy( m.value_data().begin(), m.value_data().begin() + filled2_, value_data_.begin()); 2882 index1_data_ [filled1_ - 1] = k_based(filled2_); 2883 storage_invariants (); 2884 } 2885 2886 template<class AE> 2887 BOOST_UBLAS_INLINE compressed_matrix(const matrix_expression<AE> & ae,size_type non_zeros=0)2888 compressed_matrix (const matrix_expression<AE> &ae, size_type non_zeros = 0): 2889 matrix_container<self_type> (), 2890 size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), capacity_ (restrict_capacity (non_zeros)), 2891 filled1_ (1), filled2_ (0), 2892 index1_data_ (layout_type::size_M (ae ().size1 (), ae ().size2 ()) + 1), 2893 index2_data_ (capacity_), value_data_ (capacity_) { 2894 index1_data_ [filled1_ - 1] = k_based (filled2_); 2895 storage_invariants (); 2896 matrix_assign<scalar_assign> (*this, ae); 2897 } 2898 2899 // Accessors 2900 BOOST_UBLAS_INLINE size1() const2901 size_type size1 () const { 2902 return size1_; 2903 } 2904 BOOST_UBLAS_INLINE size2() const2905 size_type size2 () const { 2906 return size2_; 2907 } 2908 BOOST_UBLAS_INLINE nnz_capacity() const2909 size_type nnz_capacity () const { 2910 return capacity_; 2911 } 2912 BOOST_UBLAS_INLINE nnz() const2913 size_type nnz () const { 2914 return filled2_; 2915 } 2916 2917 // Storage accessors 2918 BOOST_UBLAS_INLINE index_base()2919 static size_type index_base () { 2920 return IB; 2921 } 2922 BOOST_UBLAS_INLINE filled1() const2923 array_size_type filled1 () const { 2924 return filled1_; 2925 } 2926 BOOST_UBLAS_INLINE filled2() const2927 array_size_type filled2 () const { 2928 return filled2_; 2929 } 2930 BOOST_UBLAS_INLINE index1_data() const2931 const index_array_type &index1_data () const { 2932 return index1_data_; 2933 } 2934 BOOST_UBLAS_INLINE index2_data() const2935 const index_array_type &index2_data () const { 2936 return index2_data_; 2937 } 2938 BOOST_UBLAS_INLINE value_data() const2939 const value_array_type &value_data () const { 2940 return value_data_; 2941 } 2942 BOOST_UBLAS_INLINE set_filled(const array_size_type & filled1,const array_size_type & filled2)2943 void set_filled (const array_size_type& filled1, const array_size_type& filled2) { 2944 filled1_ = filled1; 2945 filled2_ = filled2; 2946 storage_invariants (); 2947 } 2948 BOOST_UBLAS_INLINE index1_data()2949 index_array_type &index1_data () { 2950 return index1_data_; 2951 } 2952 BOOST_UBLAS_INLINE index2_data()2953 index_array_type &index2_data () { 2954 return index2_data_; 2955 } 2956 BOOST_UBLAS_INLINE value_data()2957 value_array_type &value_data () { 2958 return value_data_; 2959 } 2960 BOOST_UBLAS_INLINE complete_index1_data()2961 void complete_index1_data () { 2962 while (filled1_ <= layout_type::size_M (size1_, size2_)) { 2963 this->index1_data_ [filled1_] = k_based (filled2_); 2964 ++ this->filled1_; 2965 } 2966 } 2967 2968 // Resizing 2969 private: 2970 BOOST_UBLAS_INLINE restrict_capacity(size_type non_zeros) const2971 size_type restrict_capacity (size_type non_zeros) const { 2972 non_zeros = (std::max) (non_zeros, (std::min) (size1_, size2_)); 2973 // Guarding against overflow - Thanks to Alexei Novakov for the hint. 2974 // non_zeros = (std::min) (non_zeros, size1_ * size2_); 2975 if (size1_ > 0 && non_zeros / size1_ >= size2_) 2976 non_zeros = size1_ * size2_; 2977 return non_zeros; 2978 } 2979 public: 2980 BOOST_UBLAS_INLINE resize(size_type size1,size_type size2,bool preserve=true)2981 void resize (size_type size1, size_type size2, bool preserve = true) { 2982 // FIXME preserve unimplemented 2983 BOOST_UBLAS_CHECK (!preserve, internal_logic ()); 2984 size1_ = size1; 2985 size2_ = size2; 2986 capacity_ = restrict_capacity (capacity_); 2987 filled1_ = 1; 2988 filled2_ = 0; 2989 index1_data_.resize (layout_type::size_M (size1_, size2_) + 1); 2990 index2_data_.resize (capacity_); 2991 value_data_.resize (capacity_); 2992 index1_data_ [filled1_ - 1] = k_based (filled2_); 2993 storage_invariants (); 2994 } 2995 2996 // Reserving 2997 BOOST_UBLAS_INLINE reserve(size_type non_zeros,bool preserve=true)2998 void reserve (size_type non_zeros, bool preserve = true) { 2999 capacity_ = restrict_capacity (non_zeros); 3000 if (preserve) { 3001 index2_data_.resize (capacity_, size_type ()); 3002 value_data_.resize (capacity_, value_type ()); 3003 filled2_ = (std::min) (capacity_, filled2_); 3004 } 3005 else { 3006 index2_data_.resize (capacity_); 3007 value_data_.resize (capacity_); 3008 filled1_ = 1; 3009 filled2_ = 0; 3010 index1_data_ [filled1_ - 1] = k_based (filled2_); 3011 } 3012 storage_invariants (); 3013 } 3014 3015 // Element support 3016 BOOST_UBLAS_INLINE find_element(size_type i,size_type j)3017 pointer find_element (size_type i, size_type j) { 3018 return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i, j)); 3019 } 3020 BOOST_UBLAS_INLINE find_element(size_type i,size_type j) const3021 const_pointer find_element (size_type i, size_type j) const { 3022 size_type element1 (layout_type::index_M (i, j)); 3023 size_type element2 (layout_type::index_m (i, j)); 3024 if (filled1_ <= element1 + 1) 3025 return 0; 3026 vector_const_subiterator_type itv (index1_data_.begin () + element1); 3027 const_subiterator_type it_begin (index2_data_.begin () + zero_based (*itv)); 3028 const_subiterator_type it_end (index2_data_.begin () + zero_based (*(itv + 1))); 3029 const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (element2), std::less<size_type> ())); 3030 if (it == it_end || *it != k_based (element2)) 3031 return 0; 3032 return &value_data_ [it - index2_data_.begin ()]; 3033 } 3034 3035 // Element access 3036 BOOST_UBLAS_INLINE operator ()(size_type i,size_type j) const3037 const_reference operator () (size_type i, size_type j) const { 3038 const_pointer p = find_element (i, j); 3039 if (p) 3040 return *p; 3041 else 3042 return zero_; 3043 } 3044 BOOST_UBLAS_INLINE operator ()(size_type i,size_type j)3045 reference operator () (size_type i, size_type j) { 3046 #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE 3047 size_type element1 (layout_type::index_M (i, j)); 3048 size_type element2 (layout_type::index_m (i, j)); 3049 if (filled1_ <= element1 + 1) 3050 return insert_element (i, j, value_type/*zero*/()); 3051 pointer p = find_element (i, j); 3052 if (p) 3053 return *p; 3054 else 3055 return insert_element (i, j, value_type/*zero*/()); 3056 #else 3057 return reference (*this, i, j); 3058 #endif 3059 } 3060 3061 // Element assignment 3062 BOOST_UBLAS_INLINE insert_element(size_type i,size_type j,const_reference t)3063 true_reference insert_element (size_type i, size_type j, const_reference t) { 3064 BOOST_UBLAS_CHECK (!find_element (i, j), bad_index ()); // duplicate element 3065 if (filled2_ >= capacity_) 3066 reserve (2 * filled2_, true); 3067 BOOST_UBLAS_CHECK (filled2_ < capacity_, internal_logic ()); 3068 size_type element1 = layout_type::index_M (i, j); 3069 size_type element2 = layout_type::index_m (i, j); 3070 while (filled1_ <= element1 + 1) { 3071 index1_data_ [filled1_] = k_based (filled2_); 3072 ++ filled1_; 3073 } 3074 vector_subiterator_type itv (index1_data_.begin () + element1); 3075 subiterator_type it_begin (index2_data_.begin () + zero_based (*itv)); 3076 subiterator_type it_end (index2_data_.begin () + zero_based (*(itv + 1))); 3077 subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (element2), std::less<size_type> ())); 3078 typename std::iterator_traits<subiterator_type>::difference_type n = it - index2_data_.begin (); 3079 BOOST_UBLAS_CHECK (it == it_end || *it != k_based (element2), internal_logic ()); // duplicate bound by lower_bound 3080 ++ filled2_; 3081 it = index2_data_.begin () + n; 3082 std::copy_backward (it, index2_data_.begin () + filled2_ - 1, index2_data_.begin () + filled2_); 3083 *it = k_based (element2); 3084 typename value_array_type::iterator itt (value_data_.begin () + n); 3085 std::copy_backward (itt, value_data_.begin () + filled2_ - 1, value_data_.begin () + filled2_); 3086 *itt = t; 3087 while (element1 + 1 < filled1_) { 3088 ++ index1_data_ [element1 + 1]; 3089 ++ element1; 3090 } 3091 storage_invariants (); 3092 return *itt; 3093 } 3094 BOOST_UBLAS_INLINE erase_element(size_type i,size_type j)3095 void erase_element (size_type i, size_type j) { 3096 size_type element1 = layout_type::index_M (i, j); 3097 size_type element2 = layout_type::index_m (i, j); 3098 if (element1 + 1 >= filled1_) 3099 return; 3100 vector_subiterator_type itv (index1_data_.begin () + element1); 3101 subiterator_type it_begin (index2_data_.begin () + zero_based (*itv)); 3102 subiterator_type it_end (index2_data_.begin () + zero_based (*(itv + 1))); 3103 subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (element2), std::less<size_type> ())); 3104 if (it != it_end && *it == k_based (element2)) { 3105 typename std::iterator_traits<subiterator_type>::difference_type n = it - index2_data_.begin (); 3106 std::copy (it + 1, index2_data_.begin () + filled2_, it); 3107 typename value_array_type::iterator itt (value_data_.begin () + n); 3108 std::copy (itt + 1, value_data_.begin () + filled2_, itt); 3109 -- filled2_; 3110 while (index1_data_ [filled1_ - 2] > k_based (filled2_)) { 3111 index1_data_ [filled1_ - 1] = 0; 3112 -- filled1_; 3113 } 3114 while (element1 + 1 < filled1_) { 3115 -- index1_data_ [element1 + 1]; 3116 ++ element1; 3117 } 3118 } 3119 storage_invariants (); 3120 } 3121 3122 // Zeroing 3123 BOOST_UBLAS_INLINE clear()3124 void clear () { 3125 filled1_ = 1; 3126 filled2_ = 0; 3127 index1_data_ [filled1_ - 1] = k_based (filled2_); 3128 storage_invariants (); 3129 } 3130 3131 // Assignment 3132 BOOST_UBLAS_INLINE operator =(const compressed_matrix & m)3133 compressed_matrix &operator = (const compressed_matrix &m) { 3134 if (this != &m) { 3135 size1_ = m.size1_; 3136 size2_ = m.size2_; 3137 capacity_ = m.capacity_; 3138 filled1_ = m.filled1_; 3139 filled2_ = m.filled2_; 3140 index1_data_ = m.index1_data_; 3141 index2_data_ = m.index2_data_; 3142 value_data_ = m.value_data_; 3143 } 3144 storage_invariants (); 3145 return *this; 3146 } 3147 template<class C> // Container assignment without temporary 3148 BOOST_UBLAS_INLINE operator =(const matrix_container<C> & m)3149 compressed_matrix &operator = (const matrix_container<C> &m) { 3150 resize (m ().size1 (), m ().size2 (), false); 3151 assign (m); 3152 return *this; 3153 } 3154 BOOST_UBLAS_INLINE assign_temporary(compressed_matrix & m)3155 compressed_matrix &assign_temporary (compressed_matrix &m) { 3156 swap (m); 3157 return *this; 3158 } 3159 template<class AE> 3160 BOOST_UBLAS_INLINE operator =(const matrix_expression<AE> & ae)3161 compressed_matrix &operator = (const matrix_expression<AE> &ae) { 3162 self_type temporary (ae, capacity_); 3163 return assign_temporary (temporary); 3164 } 3165 template<class AE> 3166 BOOST_UBLAS_INLINE assign(const matrix_expression<AE> & ae)3167 compressed_matrix &assign (const matrix_expression<AE> &ae) { 3168 matrix_assign<scalar_assign> (*this, ae); 3169 return *this; 3170 } 3171 template<class AE> 3172 BOOST_UBLAS_INLINE operator +=(const matrix_expression<AE> & ae)3173 compressed_matrix& operator += (const matrix_expression<AE> &ae) { 3174 self_type temporary (*this + ae, capacity_); 3175 return assign_temporary (temporary); 3176 } 3177 template<class C> // Container assignment without temporary 3178 BOOST_UBLAS_INLINE operator +=(const matrix_container<C> & m)3179 compressed_matrix &operator += (const matrix_container<C> &m) { 3180 plus_assign (m); 3181 return *this; 3182 } 3183 template<class AE> 3184 BOOST_UBLAS_INLINE plus_assign(const matrix_expression<AE> & ae)3185 compressed_matrix &plus_assign (const matrix_expression<AE> &ae) { 3186 matrix_assign<scalar_plus_assign> (*this, ae); 3187 return *this; 3188 } 3189 template<class AE> 3190 BOOST_UBLAS_INLINE operator -=(const matrix_expression<AE> & ae)3191 compressed_matrix& operator -= (const matrix_expression<AE> &ae) { 3192 self_type temporary (*this - ae, capacity_); 3193 return assign_temporary (temporary); 3194 } 3195 template<class C> // Container assignment without temporary 3196 BOOST_UBLAS_INLINE operator -=(const matrix_container<C> & m)3197 compressed_matrix &operator -= (const matrix_container<C> &m) { 3198 minus_assign (m); 3199 return *this; 3200 } 3201 template<class AE> 3202 BOOST_UBLAS_INLINE minus_assign(const matrix_expression<AE> & ae)3203 compressed_matrix &minus_assign (const matrix_expression<AE> &ae) { 3204 matrix_assign<scalar_minus_assign> (*this, ae); 3205 return *this; 3206 } 3207 template<class AT> 3208 BOOST_UBLAS_INLINE operator *=(const AT & at)3209 compressed_matrix& operator *= (const AT &at) { 3210 matrix_assign_scalar<scalar_multiplies_assign> (*this, at); 3211 return *this; 3212 } 3213 template<class AT> 3214 BOOST_UBLAS_INLINE operator /=(const AT & at)3215 compressed_matrix& operator /= (const AT &at) { 3216 matrix_assign_scalar<scalar_divides_assign> (*this, at); 3217 return *this; 3218 } 3219 3220 // Swapping 3221 BOOST_UBLAS_INLINE swap(compressed_matrix & m)3222 void swap (compressed_matrix &m) { 3223 if (this != &m) { 3224 std::swap (size1_, m.size1_); 3225 std::swap (size2_, m.size2_); 3226 std::swap (capacity_, m.capacity_); 3227 std::swap (filled1_, m.filled1_); 3228 std::swap (filled2_, m.filled2_); 3229 index1_data_.swap (m.index1_data_); 3230 index2_data_.swap (m.index2_data_); 3231 value_data_.swap (m.value_data_); 3232 } 3233 storage_invariants (); 3234 } 3235 BOOST_UBLAS_INLINE swap(compressed_matrix & m1,compressed_matrix & m2)3236 friend void swap (compressed_matrix &m1, compressed_matrix &m2) { 3237 m1.swap (m2); 3238 } 3239 3240 // Back element insertion and erasure 3241 BOOST_UBLAS_INLINE push_back(size_type i,size_type j,const_reference t)3242 void push_back (size_type i, size_type j, const_reference t) { 3243 if (filled2_ >= capacity_) 3244 reserve (2 * filled2_, true); 3245 BOOST_UBLAS_CHECK (filled2_ < capacity_, internal_logic ()); 3246 size_type element1 = layout_type::index_M (i, j); 3247 size_type element2 = layout_type::index_m (i, j); 3248 while (filled1_ < element1 + 2) { 3249 index1_data_ [filled1_] = k_based (filled2_); 3250 ++ filled1_; 3251 } 3252 // must maintain sort order 3253 BOOST_UBLAS_CHECK ((filled1_ == element1 + 2 && 3254 (filled2_ == zero_based (index1_data_ [filled1_ - 2]) || 3255 index2_data_ [filled2_ - 1] < k_based (element2))), external_logic ()); 3256 ++ filled2_; 3257 index1_data_ [filled1_ - 1] = k_based (filled2_); 3258 index2_data_ [filled2_ - 1] = k_based (element2); 3259 value_data_ [filled2_ - 1] = t; 3260 storage_invariants (); 3261 } 3262 BOOST_UBLAS_INLINE pop_back()3263 void pop_back () { 3264 BOOST_UBLAS_CHECK (filled1_ > 0 && filled2_ > 0, external_logic ()); 3265 -- filled2_; 3266 while (index1_data_ [filled1_ - 2] > k_based (filled2_)) { 3267 index1_data_ [filled1_ - 1] = 0; 3268 -- filled1_; 3269 } 3270 -- index1_data_ [filled1_ - 1]; 3271 storage_invariants (); 3272 } 3273 3274 // Iterator types 3275 private: 3276 // Use index array iterator 3277 typedef typename IA::const_iterator vector_const_subiterator_type; 3278 typedef typename IA::iterator vector_subiterator_type; 3279 typedef typename IA::const_iterator const_subiterator_type; 3280 typedef typename IA::iterator subiterator_type; 3281 3282 BOOST_UBLAS_INLINE at_element(size_type i,size_type j)3283 true_reference at_element (size_type i, size_type j) { 3284 pointer p = find_element (i, j); 3285 BOOST_UBLAS_CHECK (p, bad_index ()); 3286 return *p; 3287 } 3288 3289 public: 3290 class const_iterator1; 3291 class iterator1; 3292 class const_iterator2; 3293 class iterator2; 3294 typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1; 3295 typedef reverse_iterator_base1<iterator1> reverse_iterator1; 3296 typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2; 3297 typedef reverse_iterator_base2<iterator2> reverse_iterator2; 3298 3299 // Element lookup 3300 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. find1(int rank,size_type i,size_type j,int direction=1) const3301 const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const { 3302 for (;;) { 3303 array_size_type address1 (layout_type::index_M (i, j)); 3304 array_size_type address2 (layout_type::index_m (i, j)); 3305 vector_const_subiterator_type itv (index1_data_.begin () + (std::min) (filled1_ - 1, address1)); 3306 if (filled1_ <= address1 + 1) 3307 return const_iterator1 (*this, rank, i, j, itv, index2_data_.begin () + filled2_); 3308 3309 const_subiterator_type it_begin (index2_data_.begin () + zero_based (*itv)); 3310 const_subiterator_type it_end (index2_data_.begin () + zero_based (*(itv + 1))); 3311 3312 const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ())); 3313 if (rank == 0) 3314 return const_iterator1 (*this, rank, i, j, itv, it); 3315 if (it != it_end && zero_based (*it) == address2) 3316 return const_iterator1 (*this, rank, i, j, itv, it); 3317 if (direction > 0) { 3318 if (layout_type::fast_i ()) { 3319 if (it == it_end) 3320 return const_iterator1 (*this, rank, i, j, itv, it); 3321 i = zero_based (*it); 3322 } else { 3323 if (i >= size1_) 3324 return const_iterator1 (*this, rank, i, j, itv, it); 3325 ++ i; 3326 } 3327 } else /* if (direction < 0) */ { 3328 if (layout_type::fast_i ()) { 3329 if (it == index2_data_.begin () + zero_based (*itv)) 3330 return const_iterator1 (*this, rank, i, j, itv, it); 3331 i = zero_based (*(it - 1)); 3332 } else { 3333 if (i == 0) 3334 return const_iterator1 (*this, rank, i, j, itv, it); 3335 -- i; 3336 } 3337 } 3338 } 3339 } 3340 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. find1(int rank,size_type i,size_type j,int direction=1)3341 iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) { 3342 for (;;) { 3343 array_size_type address1 (layout_type::index_M (i, j)); 3344 array_size_type address2 (layout_type::index_m (i, j)); 3345 vector_subiterator_type itv (index1_data_.begin () + (std::min) (filled1_ - 1, address1)); 3346 if (filled1_ <= address1 + 1) 3347 return iterator1 (*this, rank, i, j, itv, index2_data_.begin () + filled2_); 3348 3349 subiterator_type it_begin (index2_data_.begin () + zero_based (*itv)); 3350 subiterator_type it_end (index2_data_.begin () + zero_based (*(itv + 1))); 3351 3352 subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ())); 3353 if (rank == 0) 3354 return iterator1 (*this, rank, i, j, itv, it); 3355 if (it != it_end && zero_based (*it) == address2) 3356 return iterator1 (*this, rank, i, j, itv, it); 3357 if (direction > 0) { 3358 if (layout_type::fast_i ()) { 3359 if (it == it_end) 3360 return iterator1 (*this, rank, i, j, itv, it); 3361 i = zero_based (*it); 3362 } else { 3363 if (i >= size1_) 3364 return iterator1 (*this, rank, i, j, itv, it); 3365 ++ i; 3366 } 3367 } else /* if (direction < 0) */ { 3368 if (layout_type::fast_i ()) { 3369 if (it == index2_data_.begin () + zero_based (*itv)) 3370 return iterator1 (*this, rank, i, j, itv, it); 3371 i = zero_based (*(it - 1)); 3372 } else { 3373 if (i == 0) 3374 return iterator1 (*this, rank, i, j, itv, it); 3375 -- i; 3376 } 3377 } 3378 } 3379 } 3380 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. find2(int rank,size_type i,size_type j,int direction=1) const3381 const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const { 3382 for (;;) { 3383 array_size_type address1 (layout_type::index_M (i, j)); 3384 array_size_type address2 (layout_type::index_m (i, j)); 3385 vector_const_subiterator_type itv (index1_data_.begin () + (std::min) (filled1_ - 1, address1)); 3386 if (filled1_ <= address1 + 1) 3387 return const_iterator2 (*this, rank, i, j, itv, index2_data_.begin () + filled2_); 3388 3389 const_subiterator_type it_begin (index2_data_.begin () + zero_based (*itv)); 3390 const_subiterator_type it_end (index2_data_.begin () + zero_based (*(itv + 1))); 3391 3392 const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ())); 3393 if (rank == 0) 3394 return const_iterator2 (*this, rank, i, j, itv, it); 3395 if (it != it_end && zero_based (*it) == address2) 3396 return const_iterator2 (*this, rank, i, j, itv, it); 3397 if (direction > 0) { 3398 if (layout_type::fast_j ()) { 3399 if (it == it_end) 3400 return const_iterator2 (*this, rank, i, j, itv, it); 3401 j = zero_based (*it); 3402 } else { 3403 if (j >= size2_) 3404 return const_iterator2 (*this, rank, i, j, itv, it); 3405 ++ j; 3406 } 3407 } else /* if (direction < 0) */ { 3408 if (layout_type::fast_j ()) { 3409 if (it == index2_data_.begin () + zero_based (*itv)) 3410 return const_iterator2 (*this, rank, i, j, itv, it); 3411 j = zero_based (*(it - 1)); 3412 } else { 3413 if (j == 0) 3414 return const_iterator2 (*this, rank, i, j, itv, it); 3415 -- j; 3416 } 3417 } 3418 } 3419 } 3420 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. find2(int rank,size_type i,size_type j,int direction=1)3421 iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) { 3422 for (;;) { 3423 array_size_type address1 (layout_type::index_M (i, j)); 3424 array_size_type address2 (layout_type::index_m (i, j)); 3425 vector_subiterator_type itv (index1_data_.begin () + (std::min) (filled1_ - 1, address1)); 3426 if (filled1_ <= address1 + 1) 3427 return iterator2 (*this, rank, i, j, itv, index2_data_.begin () + filled2_); 3428 3429 subiterator_type it_begin (index2_data_.begin () + zero_based (*itv)); 3430 subiterator_type it_end (index2_data_.begin () + zero_based (*(itv + 1))); 3431 3432 subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ())); 3433 if (rank == 0) 3434 return iterator2 (*this, rank, i, j, itv, it); 3435 if (it != it_end && zero_based (*it) == address2) 3436 return iterator2 (*this, rank, i, j, itv, it); 3437 if (direction > 0) { 3438 if (layout_type::fast_j ()) { 3439 if (it == it_end) 3440 return iterator2 (*this, rank, i, j, itv, it); 3441 j = zero_based (*it); 3442 } else { 3443 if (j >= size2_) 3444 return iterator2 (*this, rank, i, j, itv, it); 3445 ++ j; 3446 } 3447 } else /* if (direction < 0) */ { 3448 if (layout_type::fast_j ()) { 3449 if (it == index2_data_.begin () + zero_based (*itv)) 3450 return iterator2 (*this, rank, i, j, itv, it); 3451 j = zero_based (*(it - 1)); 3452 } else { 3453 if (j == 0) 3454 return iterator2 (*this, rank, i, j, itv, it); 3455 -- j; 3456 } 3457 } 3458 } 3459 } 3460 3461 3462 class const_iterator1: 3463 public container_const_reference<compressed_matrix>, 3464 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag, 3465 const_iterator1, value_type> { 3466 public: 3467 typedef typename compressed_matrix::value_type value_type; 3468 typedef typename compressed_matrix::difference_type difference_type; 3469 typedef typename compressed_matrix::const_reference reference; 3470 typedef const typename compressed_matrix::pointer pointer; 3471 3472 typedef const_iterator2 dual_iterator_type; 3473 typedef const_reverse_iterator2 dual_reverse_iterator_type; 3474 3475 // Construction and destruction 3476 BOOST_UBLAS_INLINE const_iterator1()3477 const_iterator1 (): 3478 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {} 3479 BOOST_UBLAS_INLINE const_iterator1(const self_type & m,int rank,size_type i,size_type j,const vector_const_subiterator_type & itv,const const_subiterator_type & it)3480 const_iterator1 (const self_type &m, int rank, size_type i, size_type j, const vector_const_subiterator_type &itv, const const_subiterator_type &it): 3481 container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {} 3482 BOOST_UBLAS_INLINE const_iterator1(const iterator1 & it)3483 const_iterator1 (const iterator1 &it): 3484 container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {} 3485 3486 // Arithmetic 3487 BOOST_UBLAS_INLINE operator ++()3488 const_iterator1 &operator ++ () { 3489 if (rank_ == 1 && layout_type::fast_i ()) 3490 ++ it_; 3491 else { 3492 i_ = index1 () + 1; 3493 if (rank_ == 1) 3494 *this = (*this) ().find1 (rank_, i_, j_, 1); 3495 } 3496 return *this; 3497 } 3498 BOOST_UBLAS_INLINE operator --()3499 const_iterator1 &operator -- () { 3500 if (rank_ == 1 && layout_type::fast_i ()) 3501 -- it_; 3502 else { 3503 --i_; 3504 if (rank_ == 1) 3505 *this = (*this) ().find1 (rank_, i_, j_, -1); 3506 } 3507 return *this; 3508 } 3509 3510 // Dereference 3511 BOOST_UBLAS_INLINE operator *() const3512 const_reference operator * () const { 3513 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ()); 3514 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ()); 3515 if (rank_ == 1) { 3516 return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()]; 3517 } else { 3518 return (*this) () (i_, j_); 3519 } 3520 } 3521 3522 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 3523 BOOST_UBLAS_INLINE 3524 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 3525 typename self_type:: 3526 #endif begin() const3527 const_iterator2 begin () const { 3528 const self_type &m = (*this) (); 3529 return m.find2 (1, index1 (), 0); 3530 } 3531 BOOST_UBLAS_INLINE 3532 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 3533 typename self_type:: 3534 #endif cbegin() const3535 const_iterator2 cbegin () const { 3536 return begin (); 3537 } 3538 BOOST_UBLAS_INLINE 3539 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 3540 typename self_type:: 3541 #endif end() const3542 const_iterator2 end () const { 3543 const self_type &m = (*this) (); 3544 return m.find2 (1, index1 (), m.size2 ()); 3545 } 3546 BOOST_UBLAS_INLINE 3547 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 3548 typename self_type:: 3549 #endif cend() const3550 const_iterator2 cend () const { 3551 return end (); 3552 } 3553 BOOST_UBLAS_INLINE 3554 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 3555 typename self_type:: 3556 #endif rbegin() const3557 const_reverse_iterator2 rbegin () const { 3558 return const_reverse_iterator2 (end ()); 3559 } 3560 BOOST_UBLAS_INLINE 3561 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 3562 typename self_type:: 3563 #endif crbegin() const3564 const_reverse_iterator2 crbegin () const { 3565 return rbegin (); 3566 } 3567 BOOST_UBLAS_INLINE 3568 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 3569 typename self_type:: 3570 #endif rend() const3571 const_reverse_iterator2 rend () const { 3572 return const_reverse_iterator2 (begin ()); 3573 } 3574 BOOST_UBLAS_INLINE 3575 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 3576 typename self_type:: 3577 #endif crend() const3578 const_reverse_iterator2 crend () const { 3579 return rend (); 3580 } 3581 #endif 3582 3583 // Indices 3584 BOOST_UBLAS_INLINE index1() const3585 size_type index1 () const { 3586 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ()); 3587 if (rank_ == 1) { 3588 BOOST_UBLAS_CHECK (layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ()); 3589 return layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)); 3590 } else { 3591 return i_; 3592 } 3593 } 3594 BOOST_UBLAS_INLINE index2() const3595 size_type index2 () const { 3596 if (rank_ == 1) { 3597 BOOST_UBLAS_CHECK (layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ()); 3598 return layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)); 3599 } else { 3600 return j_; 3601 } 3602 } 3603 3604 // Assignment 3605 BOOST_UBLAS_INLINE operator =(const const_iterator1 & it)3606 const_iterator1 &operator = (const const_iterator1 &it) { 3607 container_const_reference<self_type>::assign (&it ()); 3608 rank_ = it.rank_; 3609 i_ = it.i_; 3610 j_ = it.j_; 3611 itv_ = it.itv_; 3612 it_ = it.it_; 3613 return *this; 3614 } 3615 3616 // Comparison 3617 BOOST_UBLAS_INLINE operator ==(const const_iterator1 & it) const3618 bool operator == (const const_iterator1 &it) const { 3619 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 3620 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ()); 3621 if (rank_ == 1 || it.rank_ == 1) { 3622 return it_ == it.it_; 3623 } else { 3624 return i_ == it.i_ && j_ == it.j_; 3625 } 3626 } 3627 3628 private: 3629 int rank_; 3630 size_type i_; 3631 size_type j_; 3632 vector_const_subiterator_type itv_; 3633 const_subiterator_type it_; 3634 }; 3635 3636 BOOST_UBLAS_INLINE begin1() const3637 const_iterator1 begin1 () const { 3638 return find1 (0, 0, 0); 3639 } 3640 BOOST_UBLAS_INLINE cbegin1() const3641 const_iterator1 cbegin1 () const { 3642 return begin1 (); 3643 } 3644 BOOST_UBLAS_INLINE end1() const3645 const_iterator1 end1 () const { 3646 return find1 (0, size1_, 0); 3647 } 3648 BOOST_UBLAS_INLINE cend1() const3649 const_iterator1 cend1 () const { 3650 return end1 (); 3651 } 3652 3653 class iterator1: 3654 public container_reference<compressed_matrix>, 3655 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag, 3656 iterator1, value_type> { 3657 public: 3658 typedef typename compressed_matrix::value_type value_type; 3659 typedef typename compressed_matrix::difference_type difference_type; 3660 typedef typename compressed_matrix::true_reference reference; 3661 typedef typename compressed_matrix::pointer pointer; 3662 3663 typedef iterator2 dual_iterator_type; 3664 typedef reverse_iterator2 dual_reverse_iterator_type; 3665 3666 // Construction and destruction 3667 BOOST_UBLAS_INLINE iterator1()3668 iterator1 (): 3669 container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {} 3670 BOOST_UBLAS_INLINE iterator1(self_type & m,int rank,size_type i,size_type j,const vector_subiterator_type & itv,const subiterator_type & it)3671 iterator1 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it): 3672 container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {} 3673 3674 // Arithmetic 3675 BOOST_UBLAS_INLINE operator ++()3676 iterator1 &operator ++ () { 3677 if (rank_ == 1 && layout_type::fast_i ()) 3678 ++ it_; 3679 else { 3680 i_ = index1 () + 1; 3681 if (rank_ == 1) 3682 *this = (*this) ().find1 (rank_, i_, j_, 1); 3683 } 3684 return *this; 3685 } 3686 BOOST_UBLAS_INLINE operator --()3687 iterator1 &operator -- () { 3688 if (rank_ == 1 && layout_type::fast_i ()) 3689 -- it_; 3690 else { 3691 --i_; 3692 if (rank_ == 1) 3693 *this = (*this) ().find1 (rank_, i_, j_, -1); 3694 } 3695 return *this; 3696 } 3697 3698 // Dereference 3699 BOOST_UBLAS_INLINE operator *() const3700 reference operator * () const { 3701 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ()); 3702 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ()); 3703 if (rank_ == 1) { 3704 return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()]; 3705 } else { 3706 return (*this) ().at_element (i_, j_); 3707 } 3708 } 3709 3710 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 3711 BOOST_UBLAS_INLINE 3712 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 3713 typename self_type:: 3714 #endif begin() const3715 iterator2 begin () const { 3716 self_type &m = (*this) (); 3717 return m.find2 (1, index1 (), 0); 3718 } 3719 BOOST_UBLAS_INLINE 3720 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 3721 typename self_type:: 3722 #endif end() const3723 iterator2 end () const { 3724 self_type &m = (*this) (); 3725 return m.find2 (1, index1 (), m.size2 ()); 3726 } 3727 BOOST_UBLAS_INLINE 3728 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 3729 typename self_type:: 3730 #endif rbegin() const3731 reverse_iterator2 rbegin () const { 3732 return reverse_iterator2 (end ()); 3733 } 3734 BOOST_UBLAS_INLINE 3735 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 3736 typename self_type:: 3737 #endif rend() const3738 reverse_iterator2 rend () const { 3739 return reverse_iterator2 (begin ()); 3740 } 3741 #endif 3742 3743 // Indices 3744 BOOST_UBLAS_INLINE index1() const3745 size_type index1 () const { 3746 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ()); 3747 if (rank_ == 1) { 3748 BOOST_UBLAS_CHECK (layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ()); 3749 return layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)); 3750 } else { 3751 return i_; 3752 } 3753 } 3754 BOOST_UBLAS_INLINE index2() const3755 size_type index2 () const { 3756 if (rank_ == 1) { 3757 BOOST_UBLAS_CHECK (layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ()); 3758 return layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)); 3759 } else { 3760 return j_; 3761 } 3762 } 3763 3764 // Assignment 3765 BOOST_UBLAS_INLINE operator =(const iterator1 & it)3766 iterator1 &operator = (const iterator1 &it) { 3767 container_reference<self_type>::assign (&it ()); 3768 rank_ = it.rank_; 3769 i_ = it.i_; 3770 j_ = it.j_; 3771 itv_ = it.itv_; 3772 it_ = it.it_; 3773 return *this; 3774 } 3775 3776 // Comparison 3777 BOOST_UBLAS_INLINE operator ==(const iterator1 & it) const3778 bool operator == (const iterator1 &it) const { 3779 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 3780 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ()); 3781 if (rank_ == 1 || it.rank_ == 1) { 3782 return it_ == it.it_; 3783 } else { 3784 return i_ == it.i_ && j_ == it.j_; 3785 } 3786 } 3787 3788 private: 3789 int rank_; 3790 size_type i_; 3791 size_type j_; 3792 vector_subiterator_type itv_; 3793 subiterator_type it_; 3794 3795 friend class const_iterator1; 3796 }; 3797 3798 BOOST_UBLAS_INLINE begin1()3799 iterator1 begin1 () { 3800 return find1 (0, 0, 0); 3801 } 3802 BOOST_UBLAS_INLINE end1()3803 iterator1 end1 () { 3804 return find1 (0, size1_, 0); 3805 } 3806 3807 class const_iterator2: 3808 public container_const_reference<compressed_matrix>, 3809 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag, 3810 const_iterator2, value_type> { 3811 public: 3812 typedef typename compressed_matrix::value_type value_type; 3813 typedef typename compressed_matrix::difference_type difference_type; 3814 typedef typename compressed_matrix::const_reference reference; 3815 typedef const typename compressed_matrix::pointer pointer; 3816 3817 typedef const_iterator1 dual_iterator_type; 3818 typedef const_reverse_iterator1 dual_reverse_iterator_type; 3819 3820 // Construction and destruction 3821 BOOST_UBLAS_INLINE const_iterator2()3822 const_iterator2 (): 3823 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {} 3824 BOOST_UBLAS_INLINE const_iterator2(const self_type & m,int rank,size_type i,size_type j,const vector_const_subiterator_type itv,const const_subiterator_type & it)3825 const_iterator2 (const self_type &m, int rank, size_type i, size_type j, const vector_const_subiterator_type itv, const const_subiterator_type &it): 3826 container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {} 3827 BOOST_UBLAS_INLINE const_iterator2(const iterator2 & it)3828 const_iterator2 (const iterator2 &it): 3829 container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {} 3830 3831 // Arithmetic 3832 BOOST_UBLAS_INLINE operator ++()3833 const_iterator2 &operator ++ () { 3834 if (rank_ == 1 && layout_type::fast_j ()) 3835 ++ it_; 3836 else { 3837 j_ = index2 () + 1; 3838 if (rank_ == 1) 3839 *this = (*this) ().find2 (rank_, i_, j_, 1); 3840 } 3841 return *this; 3842 } 3843 BOOST_UBLAS_INLINE operator --()3844 const_iterator2 &operator -- () { 3845 if (rank_ == 1 && layout_type::fast_j ()) 3846 -- it_; 3847 else { 3848 --j_; 3849 if (rank_ == 1) 3850 *this = (*this) ().find2 (rank_, i_, j_, -1); 3851 } 3852 return *this; 3853 } 3854 3855 // Dereference 3856 BOOST_UBLAS_INLINE operator *() const3857 const_reference operator * () const { 3858 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ()); 3859 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ()); 3860 if (rank_ == 1) { 3861 return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()]; 3862 } else { 3863 return (*this) () (i_, j_); 3864 } 3865 } 3866 3867 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 3868 BOOST_UBLAS_INLINE 3869 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 3870 typename self_type:: 3871 #endif begin() const3872 const_iterator1 begin () const { 3873 const self_type &m = (*this) (); 3874 return m.find1 (1, 0, index2 ()); 3875 } 3876 BOOST_UBLAS_INLINE 3877 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 3878 typename self_type:: 3879 #endif cbegin() const3880 const_iterator1 cbegin () const { 3881 return begin (); 3882 } 3883 BOOST_UBLAS_INLINE 3884 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 3885 typename self_type:: 3886 #endif end() const3887 const_iterator1 end () const { 3888 const self_type &m = (*this) (); 3889 return m.find1 (1, m.size1 (), index2 ()); 3890 } 3891 BOOST_UBLAS_INLINE 3892 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 3893 typename self_type:: 3894 #endif cend() const3895 const_iterator1 cend () const { 3896 return end (); 3897 } 3898 BOOST_UBLAS_INLINE 3899 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 3900 typename self_type:: 3901 #endif rbegin() const3902 const_reverse_iterator1 rbegin () const { 3903 return const_reverse_iterator1 (end ()); 3904 } 3905 BOOST_UBLAS_INLINE 3906 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 3907 typename self_type:: 3908 #endif crbegin() const3909 const_reverse_iterator1 crbegin () const { 3910 return rbegin (); 3911 } 3912 BOOST_UBLAS_INLINE 3913 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 3914 typename self_type:: 3915 #endif rend() const3916 const_reverse_iterator1 rend () const { 3917 return const_reverse_iterator1 (begin ()); 3918 } 3919 BOOST_UBLAS_INLINE 3920 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 3921 typename self_type:: 3922 #endif crend() const3923 const_reverse_iterator1 crend () const { 3924 return rend (); 3925 } 3926 #endif 3927 3928 // Indices 3929 BOOST_UBLAS_INLINE index1() const3930 size_type index1 () const { 3931 if (rank_ == 1) { 3932 BOOST_UBLAS_CHECK (layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ()); 3933 return layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)); 3934 } else { 3935 return i_; 3936 } 3937 } 3938 BOOST_UBLAS_INLINE index2() const3939 size_type index2 () const { 3940 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ()); 3941 if (rank_ == 1) { 3942 BOOST_UBLAS_CHECK (layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ()); 3943 return layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)); 3944 } else { 3945 return j_; 3946 } 3947 } 3948 3949 // Assignment 3950 BOOST_UBLAS_INLINE operator =(const const_iterator2 & it)3951 const_iterator2 &operator = (const const_iterator2 &it) { 3952 container_const_reference<self_type>::assign (&it ()); 3953 rank_ = it.rank_; 3954 i_ = it.i_; 3955 j_ = it.j_; 3956 itv_ = it.itv_; 3957 it_ = it.it_; 3958 return *this; 3959 } 3960 3961 // Comparison 3962 BOOST_UBLAS_INLINE operator ==(const const_iterator2 & it) const3963 bool operator == (const const_iterator2 &it) const { 3964 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 3965 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ()); 3966 if (rank_ == 1 || it.rank_ == 1) { 3967 return it_ == it.it_; 3968 } else { 3969 return i_ == it.i_ && j_ == it.j_; 3970 } 3971 } 3972 3973 private: 3974 int rank_; 3975 size_type i_; 3976 size_type j_; 3977 vector_const_subiterator_type itv_; 3978 const_subiterator_type it_; 3979 }; 3980 3981 BOOST_UBLAS_INLINE begin2() const3982 const_iterator2 begin2 () const { 3983 return find2 (0, 0, 0); 3984 } 3985 BOOST_UBLAS_INLINE cbegin2() const3986 const_iterator2 cbegin2 () const { 3987 return begin2 (); 3988 } 3989 BOOST_UBLAS_INLINE end2() const3990 const_iterator2 end2 () const { 3991 return find2 (0, 0, size2_); 3992 } 3993 BOOST_UBLAS_INLINE cend2() const3994 const_iterator2 cend2 () const { 3995 return end2 (); 3996 } 3997 3998 class iterator2: 3999 public container_reference<compressed_matrix>, 4000 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag, 4001 iterator2, value_type> { 4002 public: 4003 typedef typename compressed_matrix::value_type value_type; 4004 typedef typename compressed_matrix::difference_type difference_type; 4005 typedef typename compressed_matrix::true_reference reference; 4006 typedef typename compressed_matrix::pointer pointer; 4007 4008 typedef iterator1 dual_iterator_type; 4009 typedef reverse_iterator1 dual_reverse_iterator_type; 4010 4011 // Construction and destruction 4012 BOOST_UBLAS_INLINE iterator2()4013 iterator2 (): 4014 container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {} 4015 BOOST_UBLAS_INLINE iterator2(self_type & m,int rank,size_type i,size_type j,const vector_subiterator_type & itv,const subiterator_type & it)4016 iterator2 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it): 4017 container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {} 4018 4019 // Arithmetic 4020 BOOST_UBLAS_INLINE operator ++()4021 iterator2 &operator ++ () { 4022 if (rank_ == 1 && layout_type::fast_j ()) 4023 ++ it_; 4024 else { 4025 j_ = index2 () + 1; 4026 if (rank_ == 1) 4027 *this = (*this) ().find2 (rank_, i_, j_, 1); 4028 } 4029 return *this; 4030 } 4031 BOOST_UBLAS_INLINE operator --()4032 iterator2 &operator -- () { 4033 if (rank_ == 1 && layout_type::fast_j ()) 4034 -- it_; 4035 else { 4036 --j_; 4037 if (rank_ == 1) 4038 *this = (*this) ().find2 (rank_, i_, j_, -1); 4039 } 4040 return *this; 4041 } 4042 4043 // Dereference 4044 BOOST_UBLAS_INLINE operator *() const4045 reference operator * () const { 4046 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ()); 4047 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ()); 4048 if (rank_ == 1) { 4049 return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()]; 4050 } else { 4051 return (*this) ().at_element (i_, j_); 4052 } 4053 } 4054 4055 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 4056 BOOST_UBLAS_INLINE 4057 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 4058 typename self_type:: 4059 #endif begin() const4060 iterator1 begin () const { 4061 self_type &m = (*this) (); 4062 return m.find1 (1, 0, index2 ()); 4063 } 4064 BOOST_UBLAS_INLINE 4065 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 4066 typename self_type:: 4067 #endif end() const4068 iterator1 end () const { 4069 self_type &m = (*this) (); 4070 return m.find1 (1, m.size1 (), index2 ()); 4071 } 4072 BOOST_UBLAS_INLINE 4073 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 4074 typename self_type:: 4075 #endif rbegin() const4076 reverse_iterator1 rbegin () const { 4077 return reverse_iterator1 (end ()); 4078 } 4079 BOOST_UBLAS_INLINE 4080 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 4081 typename self_type:: 4082 #endif rend() const4083 reverse_iterator1 rend () const { 4084 return reverse_iterator1 (begin ()); 4085 } 4086 #endif 4087 4088 // Indices 4089 BOOST_UBLAS_INLINE index1() const4090 size_type index1 () const { 4091 if (rank_ == 1) { 4092 BOOST_UBLAS_CHECK (layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ()); 4093 return layout_type::index_M (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)); 4094 } else { 4095 return i_; 4096 } 4097 } 4098 BOOST_UBLAS_INLINE index2() const4099 size_type index2 () const { 4100 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ()); 4101 if (rank_ == 1) { 4102 BOOST_UBLAS_CHECK (layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ()); 4103 return layout_type::index_m (itv_ - (*this) ().index1_data_.begin (), (*this) ().zero_based (*it_)); 4104 } else { 4105 return j_; 4106 } 4107 } 4108 4109 // Assignment 4110 BOOST_UBLAS_INLINE operator =(const iterator2 & it)4111 iterator2 &operator = (const iterator2 &it) { 4112 container_reference<self_type>::assign (&it ()); 4113 rank_ = it.rank_; 4114 i_ = it.i_; 4115 j_ = it.j_; 4116 itv_ = it.itv_; 4117 it_ = it.it_; 4118 return *this; 4119 } 4120 4121 // Comparison 4122 BOOST_UBLAS_INLINE operator ==(const iterator2 & it) const4123 bool operator == (const iterator2 &it) const { 4124 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 4125 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ()); 4126 if (rank_ == 1 || it.rank_ == 1) { 4127 return it_ == it.it_; 4128 } else { 4129 return i_ == it.i_ && j_ == it.j_; 4130 } 4131 } 4132 4133 private: 4134 int rank_; 4135 size_type i_; 4136 size_type j_; 4137 vector_subiterator_type itv_; 4138 subiterator_type it_; 4139 4140 friend class const_iterator2; 4141 }; 4142 4143 BOOST_UBLAS_INLINE begin2()4144 iterator2 begin2 () { 4145 return find2 (0, 0, 0); 4146 } 4147 BOOST_UBLAS_INLINE end2()4148 iterator2 end2 () { 4149 return find2 (0, 0, size2_); 4150 } 4151 4152 // Reverse iterators 4153 4154 BOOST_UBLAS_INLINE rbegin1() const4155 const_reverse_iterator1 rbegin1 () const { 4156 return const_reverse_iterator1 (end1 ()); 4157 } 4158 BOOST_UBLAS_INLINE crbegin1() const4159 const_reverse_iterator1 crbegin1 () const { 4160 return rbegin1 (); 4161 } 4162 BOOST_UBLAS_INLINE rend1() const4163 const_reverse_iterator1 rend1 () const { 4164 return const_reverse_iterator1 (begin1 ()); 4165 } 4166 BOOST_UBLAS_INLINE crend1() const4167 const_reverse_iterator1 crend1 () const { 4168 return rend1 (); 4169 } 4170 4171 BOOST_UBLAS_INLINE rbegin1()4172 reverse_iterator1 rbegin1 () { 4173 return reverse_iterator1 (end1 ()); 4174 } 4175 BOOST_UBLAS_INLINE rend1()4176 reverse_iterator1 rend1 () { 4177 return reverse_iterator1 (begin1 ()); 4178 } 4179 4180 BOOST_UBLAS_INLINE rbegin2() const4181 const_reverse_iterator2 rbegin2 () const { 4182 return const_reverse_iterator2 (end2 ()); 4183 } 4184 BOOST_UBLAS_INLINE crbegin2() const4185 const_reverse_iterator2 crbegin2 () const { 4186 return rbegin2 (); 4187 } 4188 BOOST_UBLAS_INLINE rend2() const4189 const_reverse_iterator2 rend2 () const { 4190 return const_reverse_iterator2 (begin2 ()); 4191 } 4192 BOOST_UBLAS_INLINE crend2() const4193 const_reverse_iterator2 crend2 () const { 4194 return rend2 (); 4195 } 4196 4197 BOOST_UBLAS_INLINE rbegin2()4198 reverse_iterator2 rbegin2 () { 4199 return reverse_iterator2 (end2 ()); 4200 } 4201 BOOST_UBLAS_INLINE rend2()4202 reverse_iterator2 rend2 () { 4203 return reverse_iterator2 (begin2 ()); 4204 } 4205 4206 // Serialization 4207 template<class Archive> serialize(Archive & ar,const unsigned int)4208 void serialize(Archive & ar, const unsigned int /* file_version */){ 4209 serialization::collection_size_type s1 (size1_); 4210 serialization::collection_size_type s2 (size2_); 4211 ar & serialization::make_nvp("size1",s1); 4212 ar & serialization::make_nvp("size2",s2); 4213 if (Archive::is_loading::value) { 4214 size1_ = s1; 4215 size2_ = s2; 4216 } 4217 ar & serialization::make_nvp("capacity", capacity_); 4218 ar & serialization::make_nvp("filled1", filled1_); 4219 ar & serialization::make_nvp("filled2", filled2_); 4220 ar & serialization::make_nvp("index1_data", index1_data_); 4221 ar & serialization::make_nvp("index2_data", index2_data_); 4222 ar & serialization::make_nvp("value_data", value_data_); 4223 storage_invariants(); 4224 } 4225 4226 private: storage_invariants() const4227 void storage_invariants () const { 4228 BOOST_UBLAS_CHECK (layout_type::size_M (size1_, size2_) + 1 == index1_data_.size (), internal_logic ()); 4229 BOOST_UBLAS_CHECK (capacity_ == index2_data_.size (), internal_logic ()); 4230 BOOST_UBLAS_CHECK (capacity_ == value_data_.size (), internal_logic ()); 4231 BOOST_UBLAS_CHECK (filled1_ > 0 && filled1_ <= layout_type::size_M (size1_, size2_) + 1, internal_logic ()); 4232 BOOST_UBLAS_CHECK (filled2_ <= capacity_, internal_logic ()); 4233 BOOST_UBLAS_CHECK (index1_data_ [filled1_ - 1] == k_based (filled2_), internal_logic ()); 4234 } 4235 4236 size_type size1_; 4237 size_type size2_; 4238 array_size_type capacity_; 4239 array_size_type filled1_; 4240 array_size_type filled2_; 4241 index_array_type index1_data_; 4242 index_array_type index2_data_; 4243 value_array_type value_data_; 4244 static const value_type zero_; 4245 4246 BOOST_UBLAS_INLINE zero_based(size_type k_based_index)4247 static size_type zero_based (size_type k_based_index) { 4248 return k_based_index - IB; 4249 } 4250 BOOST_UBLAS_INLINE k_based(size_type zero_based_index)4251 static size_type k_based (size_type zero_based_index) { 4252 return zero_based_index + IB; 4253 } 4254 4255 friend class iterator1; 4256 friend class iterator2; 4257 friend class const_iterator1; 4258 friend class const_iterator2; 4259 }; 4260 4261 template<class T, class L, std::size_t IB, class IA, class TA> 4262 const typename compressed_matrix<T, L, IB, IA, TA>::value_type compressed_matrix<T, L, IB, IA, TA>::zero_ = value_type/*zero*/(); 4263 4264 4265 // Coordinate array based sparse matrix class 4266 // Thanks to Kresimir Fresl for extending this to cover different index bases. 4267 template<class T, class L, std::size_t IB, class IA, class TA> 4268 class coordinate_matrix: 4269 public matrix_container<coordinate_matrix<T, L, IB, IA, TA> > { 4270 4271 typedef T &true_reference; 4272 typedef T *pointer; 4273 typedef const T *const_pointer; 4274 typedef L layout_type; 4275 typedef coordinate_matrix<T, L, IB, IA, TA> self_type; 4276 public: 4277 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS 4278 using matrix_container<self_type>::operator (); 4279 #endif 4280 // ISSUE require type consistency check, is_convertable (IA::size_type, TA::size_type) 4281 typedef typename IA::value_type size_type; 4282 // ISSUE difference_type cannot be deduced for sparse indices, we only know the value_type 4283 typedef std::ptrdiff_t difference_type; 4284 // size_type for the data arrays. 4285 typedef typename IA::size_type array_size_type; 4286 typedef T value_type; 4287 typedef const T &const_reference; 4288 #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE 4289 typedef T &reference; 4290 #else 4291 typedef sparse_matrix_element<self_type> reference; 4292 #endif 4293 typedef IA index_array_type; 4294 typedef TA value_array_type; 4295 typedef const matrix_reference<const self_type> const_closure_type; 4296 typedef matrix_reference<self_type> closure_type; 4297 typedef coordinate_vector<T, IB, IA, TA> vector_temporary_type; 4298 typedef self_type matrix_temporary_type; 4299 typedef sparse_tag storage_category; 4300 typedef typename L::orientation_category orientation_category; 4301 4302 // Construction and destruction 4303 BOOST_UBLAS_INLINE coordinate_matrix()4304 coordinate_matrix (): 4305 matrix_container<self_type> (), 4306 size1_ (0), size2_ (0), capacity_ (restrict_capacity (0)), 4307 filled_ (0), sorted_filled_ (filled_), sorted_ (true), 4308 index1_data_ (capacity_), index2_data_ (capacity_), value_data_ (capacity_) { 4309 storage_invariants (); 4310 } 4311 BOOST_UBLAS_INLINE coordinate_matrix(size_type size1,size_type size2,array_size_type non_zeros=0)4312 coordinate_matrix (size_type size1, size_type size2, array_size_type non_zeros = 0): 4313 matrix_container<self_type> (), 4314 size1_ (size1), size2_ (size2), capacity_ (restrict_capacity (non_zeros)), 4315 filled_ (0), sorted_filled_ (filled_), sorted_ (true), 4316 index1_data_ (capacity_), index2_data_ (capacity_), value_data_ (capacity_) { 4317 storage_invariants (); 4318 } 4319 BOOST_UBLAS_INLINE coordinate_matrix(const coordinate_matrix & m)4320 coordinate_matrix (const coordinate_matrix &m): 4321 matrix_container<self_type> (), 4322 size1_ (m.size1_), size2_ (m.size2_), capacity_ (m.capacity_), 4323 filled_ (m.filled_), sorted_filled_ (m.sorted_filled_), sorted_ (m.sorted_), 4324 index1_data_ (m.index1_data_), index2_data_ (m.index2_data_), value_data_ (m.value_data_) { 4325 storage_invariants (); 4326 } 4327 template<class AE> 4328 BOOST_UBLAS_INLINE coordinate_matrix(const matrix_expression<AE> & ae,array_size_type non_zeros=0)4329 coordinate_matrix (const matrix_expression<AE> &ae, array_size_type non_zeros = 0): 4330 matrix_container<self_type> (), 4331 size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), capacity_ (restrict_capacity (non_zeros)), 4332 filled_ (0), sorted_filled_ (filled_), sorted_ (true), 4333 index1_data_ (capacity_), index2_data_ (capacity_), value_data_ (capacity_) { 4334 storage_invariants (); 4335 matrix_assign<scalar_assign> (*this, ae); 4336 } 4337 4338 // Accessors 4339 BOOST_UBLAS_INLINE size1() const4340 size_type size1 () const { 4341 return size1_; 4342 } 4343 BOOST_UBLAS_INLINE size2() const4344 size_type size2 () const { 4345 return size2_; 4346 } 4347 BOOST_UBLAS_INLINE nnz_capacity() const4348 size_type nnz_capacity () const { 4349 return capacity_; 4350 } 4351 BOOST_UBLAS_INLINE nnz() const4352 size_type nnz () const { 4353 return filled_; 4354 } 4355 4356 // Storage accessors 4357 BOOST_UBLAS_INLINE index_base()4358 static size_type index_base () { 4359 return IB; 4360 } 4361 BOOST_UBLAS_INLINE filled() const4362 array_size_type filled () const { 4363 return filled_; 4364 } 4365 BOOST_UBLAS_INLINE index1_data() const4366 const index_array_type &index1_data () const { 4367 return index1_data_; 4368 } 4369 BOOST_UBLAS_INLINE index2_data() const4370 const index_array_type &index2_data () const { 4371 return index2_data_; 4372 } 4373 BOOST_UBLAS_INLINE value_data() const4374 const value_array_type &value_data () const { 4375 return value_data_; 4376 } 4377 BOOST_UBLAS_INLINE set_filled(const array_size_type & filled)4378 void set_filled (const array_size_type &filled) { 4379 // Make sure that storage_invariants() succeeds 4380 if (sorted_ && filled < filled_) 4381 sorted_filled_ = filled; 4382 else 4383 sorted_ = (sorted_filled_ == filled); 4384 filled_ = filled; 4385 storage_invariants (); 4386 } 4387 BOOST_UBLAS_INLINE index1_data()4388 index_array_type &index1_data () { 4389 return index1_data_; 4390 } 4391 BOOST_UBLAS_INLINE index2_data()4392 index_array_type &index2_data () { 4393 return index2_data_; 4394 } 4395 BOOST_UBLAS_INLINE value_data()4396 value_array_type &value_data () { 4397 return value_data_; 4398 } 4399 4400 // Resizing 4401 private: 4402 BOOST_UBLAS_INLINE restrict_capacity(array_size_type non_zeros) const4403 array_size_type restrict_capacity (array_size_type non_zeros) const { 4404 // minimum non_zeros 4405 non_zeros = (std::max) (non_zeros, array_size_type((std::min) (size1_, size2_))); 4406 // ISSUE no maximum as coordinate may contain inserted duplicates 4407 return non_zeros; 4408 } 4409 public: 4410 BOOST_UBLAS_INLINE resize(size_type size1,size_type size2,bool preserve=true)4411 void resize (size_type size1, size_type size2, bool preserve = true) { 4412 // FIXME preserve unimplemented 4413 BOOST_UBLAS_CHECK (!preserve, internal_logic ()); 4414 size1_ = size1; 4415 size2_ = size2; 4416 capacity_ = restrict_capacity (capacity_); 4417 index1_data_.resize (capacity_); 4418 index2_data_.resize (capacity_); 4419 value_data_.resize (capacity_); 4420 filled_ = 0; 4421 sorted_filled_ = filled_; 4422 sorted_ = true; 4423 storage_invariants (); 4424 } 4425 4426 // Reserving 4427 BOOST_UBLAS_INLINE reserve(array_size_type non_zeros,bool preserve=true)4428 void reserve (array_size_type non_zeros, bool preserve = true) { 4429 sort (); // remove duplicate elements 4430 capacity_ = restrict_capacity (non_zeros); 4431 if (preserve) { 4432 index1_data_.resize (capacity_, size_type ()); 4433 index2_data_.resize (capacity_, size_type ()); 4434 value_data_.resize (capacity_, value_type ()); 4435 filled_ = (std::min) (capacity_, filled_); 4436 } 4437 else { 4438 index1_data_.resize (capacity_); 4439 index2_data_.resize (capacity_); 4440 value_data_.resize (capacity_); 4441 filled_ = 0; 4442 } 4443 sorted_filled_ = filled_; 4444 storage_invariants (); 4445 } 4446 4447 // Element support 4448 BOOST_UBLAS_INLINE find_element(size_type i,size_type j)4449 pointer find_element (size_type i, size_type j) { 4450 return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i, j)); 4451 } 4452 BOOST_UBLAS_INLINE find_element(size_type i,size_type j) const4453 const_pointer find_element (size_type i, size_type j) const { 4454 sort (); 4455 size_type element1 (layout_type::index_M (i, j)); 4456 size_type element2 (layout_type::index_m (i, j)); 4457 vector_const_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (element1), std::less<size_type> ())); 4458 vector_const_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (element1), std::less<size_type> ())); 4459 if (itv_begin == itv_end) 4460 return 0; 4461 const_subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ())); 4462 const_subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ())); 4463 const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (element2), std::less<size_type> ())); 4464 if (it == it_end || *it != k_based (element2)) 4465 return 0; 4466 return &value_data_ [it - index2_data_.begin ()]; 4467 } 4468 4469 // Element access 4470 BOOST_UBLAS_INLINE operator ()(size_type i,size_type j) const4471 const_reference operator () (size_type i, size_type j) const { 4472 const_pointer p = find_element (i, j); 4473 if (p) 4474 return *p; 4475 else 4476 return zero_; 4477 } 4478 BOOST_UBLAS_INLINE operator ()(size_type i,size_type j)4479 reference operator () (size_type i, size_type j) { 4480 #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE 4481 pointer p = find_element (i, j); 4482 if (p) 4483 return *p; 4484 else 4485 return insert_element (i, j, value_type/*zero*/()); 4486 #else 4487 return reference (*this, i, j); 4488 #endif 4489 } 4490 4491 // Element assignment 4492 BOOST_UBLAS_INLINE append_element(size_type i,size_type j,const_reference t)4493 void append_element (size_type i, size_type j, const_reference t) { 4494 if (filled_ >= capacity_) 4495 reserve (2 * filled_, true); 4496 BOOST_UBLAS_CHECK (filled_ < capacity_, internal_logic ()); 4497 size_type element1 = layout_type::index_M (i, j); 4498 size_type element2 = layout_type::index_m (i, j); 4499 index1_data_ [filled_] = k_based (element1); 4500 index2_data_ [filled_] = k_based (element2); 4501 value_data_ [filled_] = t; 4502 ++ filled_; 4503 sorted_ = false; 4504 storage_invariants (); 4505 } 4506 BOOST_UBLAS_INLINE insert_element(size_type i,size_type j,const_reference t)4507 true_reference insert_element (size_type i, size_type j, const_reference t) { 4508 BOOST_UBLAS_CHECK (!find_element (i, j), bad_index ()); // duplicate element 4509 append_element (i, j, t); 4510 return value_data_ [filled_ - 1]; 4511 } 4512 BOOST_UBLAS_INLINE erase_element(size_type i,size_type j)4513 void erase_element (size_type i, size_type j) { 4514 size_type element1 = layout_type::index_M (i, j); 4515 size_type element2 = layout_type::index_m (i, j); 4516 sort (); 4517 vector_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (element1), std::less<size_type> ())); 4518 vector_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (element1), std::less<size_type> ())); 4519 subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ())); 4520 subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ())); 4521 subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (element2), std::less<size_type> ())); 4522 if (it != it_end && *it == k_based (element2)) { 4523 typename std::iterator_traits<subiterator_type>::difference_type n = it - index2_data_.begin (); 4524 vector_subiterator_type itv (index1_data_.begin () + n); 4525 std::copy (itv + 1, index1_data_.begin () + filled_, itv); 4526 std::copy (it + 1, index2_data_.begin () + filled_, it); 4527 typename value_array_type::iterator itt (value_data_.begin () + n); 4528 std::copy (itt + 1, value_data_.begin () + filled_, itt); 4529 -- filled_; 4530 sorted_filled_ = filled_; 4531 } 4532 storage_invariants (); 4533 } 4534 4535 // Zeroing 4536 BOOST_UBLAS_INLINE clear()4537 void clear () { 4538 filled_ = 0; 4539 sorted_filled_ = filled_; 4540 sorted_ = true; 4541 storage_invariants (); 4542 } 4543 4544 // Assignment 4545 BOOST_UBLAS_INLINE operator =(const coordinate_matrix & m)4546 coordinate_matrix &operator = (const coordinate_matrix &m) { 4547 if (this != &m) { 4548 size1_ = m.size1_; 4549 size2_ = m.size2_; 4550 capacity_ = m.capacity_; 4551 filled_ = m.filled_; 4552 sorted_filled_ = m.sorted_filled_; 4553 sorted_ = m.sorted_; 4554 index1_data_ = m.index1_data_; 4555 index2_data_ = m.index2_data_; 4556 value_data_ = m.value_data_; 4557 BOOST_UBLAS_CHECK (capacity_ == index1_data_.size (), internal_logic ()); 4558 BOOST_UBLAS_CHECK (capacity_ == index2_data_.size (), internal_logic ()); 4559 BOOST_UBLAS_CHECK (capacity_ == value_data_.size (), internal_logic ()); 4560 } 4561 storage_invariants (); 4562 return *this; 4563 } 4564 template<class C> // Container assignment without temporary 4565 BOOST_UBLAS_INLINE operator =(const matrix_container<C> & m)4566 coordinate_matrix &operator = (const matrix_container<C> &m) { 4567 resize (m ().size1 (), m ().size2 (), false); 4568 assign (m); 4569 return *this; 4570 } 4571 BOOST_UBLAS_INLINE assign_temporary(coordinate_matrix & m)4572 coordinate_matrix &assign_temporary (coordinate_matrix &m) { 4573 swap (m); 4574 return *this; 4575 } 4576 template<class AE> 4577 BOOST_UBLAS_INLINE operator =(const matrix_expression<AE> & ae)4578 coordinate_matrix &operator = (const matrix_expression<AE> &ae) { 4579 self_type temporary (ae, capacity_); 4580 return assign_temporary (temporary); 4581 } 4582 template<class AE> 4583 BOOST_UBLAS_INLINE assign(const matrix_expression<AE> & ae)4584 coordinate_matrix &assign (const matrix_expression<AE> &ae) { 4585 matrix_assign<scalar_assign> (*this, ae); 4586 return *this; 4587 } 4588 template<class AE> 4589 BOOST_UBLAS_INLINE operator +=(const matrix_expression<AE> & ae)4590 coordinate_matrix& operator += (const matrix_expression<AE> &ae) { 4591 self_type temporary (*this + ae, capacity_); 4592 return assign_temporary (temporary); 4593 } 4594 template<class C> // Container assignment without temporary 4595 BOOST_UBLAS_INLINE operator +=(const matrix_container<C> & m)4596 coordinate_matrix &operator += (const matrix_container<C> &m) { 4597 plus_assign (m); 4598 return *this; 4599 } 4600 template<class AE> 4601 BOOST_UBLAS_INLINE plus_assign(const matrix_expression<AE> & ae)4602 coordinate_matrix &plus_assign (const matrix_expression<AE> &ae) { 4603 matrix_assign<scalar_plus_assign> (*this, ae); 4604 return *this; 4605 } 4606 template<class AE> 4607 BOOST_UBLAS_INLINE operator -=(const matrix_expression<AE> & ae)4608 coordinate_matrix& operator -= (const matrix_expression<AE> &ae) { 4609 self_type temporary (*this - ae, capacity_); 4610 return assign_temporary (temporary); 4611 } 4612 template<class C> // Container assignment without temporary 4613 BOOST_UBLAS_INLINE operator -=(const matrix_container<C> & m)4614 coordinate_matrix &operator -= (const matrix_container<C> &m) { 4615 minus_assign (m); 4616 return *this; 4617 } 4618 template<class AE> 4619 BOOST_UBLAS_INLINE minus_assign(const matrix_expression<AE> & ae)4620 coordinate_matrix &minus_assign (const matrix_expression<AE> &ae) { 4621 matrix_assign<scalar_minus_assign> (*this, ae); 4622 return *this; 4623 } 4624 template<class AT> 4625 BOOST_UBLAS_INLINE operator *=(const AT & at)4626 coordinate_matrix& operator *= (const AT &at) { 4627 matrix_assign_scalar<scalar_multiplies_assign> (*this, at); 4628 return *this; 4629 } 4630 template<class AT> 4631 BOOST_UBLAS_INLINE operator /=(const AT & at)4632 coordinate_matrix& operator /= (const AT &at) { 4633 matrix_assign_scalar<scalar_divides_assign> (*this, at); 4634 return *this; 4635 } 4636 4637 // Swapping 4638 BOOST_UBLAS_INLINE swap(coordinate_matrix & m)4639 void swap (coordinate_matrix &m) { 4640 if (this != &m) { 4641 std::swap (size1_, m.size1_); 4642 std::swap (size2_, m.size2_); 4643 std::swap (capacity_, m.capacity_); 4644 std::swap (filled_, m.filled_); 4645 std::swap (sorted_filled_, m.sorted_filled_); 4646 std::swap (sorted_, m.sorted_); 4647 index1_data_.swap (m.index1_data_); 4648 index2_data_.swap (m.index2_data_); 4649 value_data_.swap (m.value_data_); 4650 } 4651 storage_invariants (); 4652 } 4653 BOOST_UBLAS_INLINE swap(coordinate_matrix & m1,coordinate_matrix & m2)4654 friend void swap (coordinate_matrix &m1, coordinate_matrix &m2) { 4655 m1.swap (m2); 4656 } 4657 4658 // replacement if STL lower bound algorithm for use of inplace_merge lower_bound(array_size_type beg,array_size_type end,array_size_type target) const4659 array_size_type lower_bound (array_size_type beg, array_size_type end, array_size_type target) const { 4660 while (end > beg) { 4661 array_size_type mid = (beg + end) / 2; 4662 if (((index1_data_[mid] < index1_data_[target]) || 4663 ((index1_data_[mid] == index1_data_[target]) && 4664 (index2_data_[mid] < index2_data_[target])))) { 4665 beg = mid + 1; 4666 } else { 4667 end = mid; 4668 } 4669 } 4670 return beg; 4671 } 4672 4673 // specialized replacement of STL inplace_merge to avoid compilation 4674 // problems with respect to the array_triple iterator inplace_merge(array_size_type beg,array_size_type mid,array_size_type end) const4675 void inplace_merge (array_size_type beg, array_size_type mid, array_size_type end) const { 4676 array_size_type len_lef = mid - beg; 4677 array_size_type len_rig = end - mid; 4678 4679 if (len_lef == 1 && len_rig == 1) { 4680 if ((index1_data_[mid] < index1_data_[beg]) || 4681 ((index1_data_[mid] == index1_data_[beg]) && (index2_data_[mid] < index2_data_[beg]))) 4682 { 4683 std::swap(index1_data_[beg], index1_data_[mid]); 4684 std::swap(index2_data_[beg], index2_data_[mid]); 4685 std::swap(value_data_[beg], value_data_[mid]); 4686 } 4687 } else if (len_lef > 0 && len_rig > 0) { 4688 array_size_type lef_mid, rig_mid; 4689 if (len_lef >= len_rig) { 4690 lef_mid = (beg + mid) / 2; 4691 rig_mid = lower_bound(mid, end, lef_mid); 4692 } else { 4693 rig_mid = (mid + end) / 2; 4694 lef_mid = lower_bound(beg, mid, rig_mid); 4695 } 4696 std::rotate(&index1_data_[0] + lef_mid, &index1_data_[0] + mid, &index1_data_[0] + rig_mid); 4697 std::rotate(&index2_data_[0] + lef_mid, &index2_data_[0] + mid, &index2_data_[0] + rig_mid); 4698 std::rotate(&value_data_[0] + lef_mid, &value_data_[0] + mid, &value_data_[0] + rig_mid); 4699 4700 array_size_type new_mid = lef_mid + rig_mid - mid; 4701 inplace_merge(beg, lef_mid, new_mid); 4702 inplace_merge(new_mid, rig_mid, end); 4703 } 4704 } 4705 4706 // Sorting and summation of duplicates 4707 BOOST_UBLAS_INLINE sort() const4708 void sort () const { 4709 if (! sorted_ && filled_ > 0) { 4710 typedef index_triple_array<index_array_type, index_array_type, value_array_type> array_triple; 4711 array_triple ita (filled_, index1_data_, index2_data_, value_data_); 4712 #ifndef BOOST_UBLAS_COO_ALWAYS_DO_FULL_SORT 4713 const typename array_triple::iterator iunsorted = ita.begin () + sorted_filled_; 4714 // sort new elements and merge 4715 std::sort (iunsorted, ita.end ()); 4716 inplace_merge(0, sorted_filled_, filled_); 4717 #else 4718 const typename array_triple::iterator iunsorted = ita.begin (); 4719 std::sort (iunsorted, ita.end ()); 4720 #endif 4721 // sum duplicates with += and remove 4722 array_size_type filled = 0; 4723 for (array_size_type i = 1; i < filled_; ++ i) { 4724 if (index1_data_ [filled] != index1_data_ [i] || 4725 index2_data_ [filled] != index2_data_ [i]) { 4726 ++ filled; 4727 if (filled != i) { 4728 index1_data_ [filled] = index1_data_ [i]; 4729 index2_data_ [filled] = index2_data_ [i]; 4730 value_data_ [filled] = value_data_ [i]; 4731 } 4732 } else { 4733 value_data_ [filled] += value_data_ [i]; 4734 } 4735 } 4736 filled_ = filled + 1; 4737 sorted_filled_ = filled_; 4738 sorted_ = true; 4739 storage_invariants (); 4740 } 4741 } 4742 4743 // Back element insertion and erasure 4744 BOOST_UBLAS_INLINE push_back(size_type i,size_type j,const_reference t)4745 void push_back (size_type i, size_type j, const_reference t) { 4746 size_type element1 = layout_type::index_M (i, j); 4747 size_type element2 = layout_type::index_m (i, j); 4748 // must maintain sort order 4749 BOOST_UBLAS_CHECK (sorted_ && 4750 (filled_ == 0 || 4751 index1_data_ [filled_ - 1] < k_based (element1) || 4752 (index1_data_ [filled_ - 1] == k_based (element1) && index2_data_ [filled_ - 1] < k_based (element2))) 4753 , external_logic ()); 4754 if (filled_ >= capacity_) 4755 reserve (2 * filled_, true); 4756 BOOST_UBLAS_CHECK (filled_ < capacity_, internal_logic ()); 4757 index1_data_ [filled_] = k_based (element1); 4758 index2_data_ [filled_] = k_based (element2); 4759 value_data_ [filled_] = t; 4760 ++ filled_; 4761 sorted_filled_ = filled_; 4762 storage_invariants (); 4763 } 4764 BOOST_UBLAS_INLINE pop_back()4765 void pop_back () { 4766 // ISSUE invariants could be simpilfied if sorted required as precondition 4767 BOOST_UBLAS_CHECK (filled_ > 0, external_logic ()); 4768 -- filled_; 4769 sorted_filled_ = (std::min) (sorted_filled_, filled_); 4770 sorted_ = sorted_filled_ = filled_; 4771 storage_invariants (); 4772 } 4773 4774 // Iterator types 4775 private: 4776 // Use index array iterator 4777 typedef typename IA::const_iterator vector_const_subiterator_type; 4778 typedef typename IA::iterator vector_subiterator_type; 4779 typedef typename IA::const_iterator const_subiterator_type; 4780 typedef typename IA::iterator subiterator_type; 4781 4782 BOOST_UBLAS_INLINE at_element(size_type i,size_type j)4783 true_reference at_element (size_type i, size_type j) { 4784 pointer p = find_element (i, j); 4785 BOOST_UBLAS_CHECK (p, bad_index ()); 4786 return *p; 4787 } 4788 4789 public: 4790 class const_iterator1; 4791 class iterator1; 4792 class const_iterator2; 4793 class iterator2; 4794 typedef reverse_iterator_base1<const_iterator1> const_reverse_iterator1; 4795 typedef reverse_iterator_base1<iterator1> reverse_iterator1; 4796 typedef reverse_iterator_base2<const_iterator2> const_reverse_iterator2; 4797 typedef reverse_iterator_base2<iterator2> reverse_iterator2; 4798 4799 // Element lookup 4800 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. find1(int rank,size_type i,size_type j,int direction=1) const4801 const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const { 4802 sort (); 4803 for (;;) { 4804 size_type address1 (layout_type::index_M (i, j)); 4805 size_type address2 (layout_type::index_m (i, j)); 4806 vector_const_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ())); 4807 vector_const_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ())); 4808 4809 const_subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ())); 4810 const_subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ())); 4811 4812 const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ())); 4813 vector_const_subiterator_type itv (index1_data_.begin () + (it - index2_data_.begin ())); 4814 if (rank == 0) 4815 return const_iterator1 (*this, rank, i, j, itv, it); 4816 if (it != it_end && zero_based (*it) == address2) 4817 return const_iterator1 (*this, rank, i, j, itv, it); 4818 if (direction > 0) { 4819 if (layout_type::fast_i ()) { 4820 if (it == it_end) 4821 return const_iterator1 (*this, rank, i, j, itv, it); 4822 i = zero_based (*it); 4823 } else { 4824 if (i >= size1_) 4825 return const_iterator1 (*this, rank, i, j, itv, it); 4826 ++ i; 4827 } 4828 } else /* if (direction < 0) */ { 4829 if (layout_type::fast_i ()) { 4830 if (it == index2_data_.begin () + array_size_type (zero_based (*itv))) 4831 return const_iterator1 (*this, rank, i, j, itv, it); 4832 i = zero_based (*(it - 1)); 4833 } else { 4834 if (i == 0) 4835 return const_iterator1 (*this, rank, i, j, itv, it); 4836 -- i; 4837 } 4838 } 4839 } 4840 } 4841 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. find1(int rank,size_type i,size_type j,int direction=1)4842 iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) { 4843 sort (); 4844 for (;;) { 4845 size_type address1 (layout_type::index_M (i, j)); 4846 size_type address2 (layout_type::index_m (i, j)); 4847 vector_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ())); 4848 vector_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ())); 4849 4850 subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ())); 4851 subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ())); 4852 4853 subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ())); 4854 vector_subiterator_type itv (index1_data_.begin () + (it - index2_data_.begin ())); 4855 if (rank == 0) 4856 return iterator1 (*this, rank, i, j, itv, it); 4857 if (it != it_end && zero_based (*it) == address2) 4858 return iterator1 (*this, rank, i, j, itv, it); 4859 if (direction > 0) { 4860 if (layout_type::fast_i ()) { 4861 if (it == it_end) 4862 return iterator1 (*this, rank, i, j, itv, it); 4863 i = zero_based (*it); 4864 } else { 4865 if (i >= size1_) 4866 return iterator1 (*this, rank, i, j, itv, it); 4867 ++ i; 4868 } 4869 } else /* if (direction < 0) */ { 4870 if (layout_type::fast_i ()) { 4871 if (it == index2_data_.begin () + array_size_type (zero_based (*itv))) 4872 return iterator1 (*this, rank, i, j, itv, it); 4873 i = zero_based (*(it - 1)); 4874 } else { 4875 if (i == 0) 4876 return iterator1 (*this, rank, i, j, itv, it); 4877 -- i; 4878 } 4879 } 4880 } 4881 } 4882 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. find2(int rank,size_type i,size_type j,int direction=1) const4883 const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const { 4884 sort (); 4885 for (;;) { 4886 size_type address1 (layout_type::index_M (i, j)); 4887 size_type address2 (layout_type::index_m (i, j)); 4888 vector_const_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ())); 4889 vector_const_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ())); 4890 4891 const_subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ())); 4892 const_subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ())); 4893 4894 const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ())); 4895 vector_const_subiterator_type itv (index1_data_.begin () + (it - index2_data_.begin ())); 4896 if (rank == 0) 4897 return const_iterator2 (*this, rank, i, j, itv, it); 4898 if (it != it_end && zero_based (*it) == address2) 4899 return const_iterator2 (*this, rank, i, j, itv, it); 4900 if (direction > 0) { 4901 if (layout_type::fast_j ()) { 4902 if (it == it_end) 4903 return const_iterator2 (*this, rank, i, j, itv, it); 4904 j = zero_based (*it); 4905 } else { 4906 if (j >= size2_) 4907 return const_iterator2 (*this, rank, i, j, itv, it); 4908 ++ j; 4909 } 4910 } else /* if (direction < 0) */ { 4911 if (layout_type::fast_j ()) { 4912 if (it == index2_data_.begin () + array_size_type (zero_based (*itv))) 4913 return const_iterator2 (*this, rank, i, j, itv, it); 4914 j = zero_based (*(it - 1)); 4915 } else { 4916 if (j == 0) 4917 return const_iterator2 (*this, rank, i, j, itv, it); 4918 -- j; 4919 } 4920 } 4921 } 4922 } 4923 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. find2(int rank,size_type i,size_type j,int direction=1)4924 iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) { 4925 sort (); 4926 for (;;) { 4927 size_type address1 (layout_type::index_M (i, j)); 4928 size_type address2 (layout_type::index_m (i, j)); 4929 vector_subiterator_type itv_begin (detail::lower_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ())); 4930 vector_subiterator_type itv_end (detail::upper_bound (index1_data_.begin (), index1_data_.begin () + filled_, k_based (address1), std::less<size_type> ())); 4931 4932 subiterator_type it_begin (index2_data_.begin () + (itv_begin - index1_data_.begin ())); 4933 subiterator_type it_end (index2_data_.begin () + (itv_end - index1_data_.begin ())); 4934 4935 subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ())); 4936 vector_subiterator_type itv (index1_data_.begin () + (it - index2_data_.begin ())); 4937 if (rank == 0) 4938 return iterator2 (*this, rank, i, j, itv, it); 4939 if (it != it_end && zero_based (*it) == address2) 4940 return iterator2 (*this, rank, i, j, itv, it); 4941 if (direction > 0) { 4942 if (layout_type::fast_j ()) { 4943 if (it == it_end) 4944 return iterator2 (*this, rank, i, j, itv, it); 4945 j = zero_based (*it); 4946 } else { 4947 if (j >= size2_) 4948 return iterator2 (*this, rank, i, j, itv, it); 4949 ++ j; 4950 } 4951 } else /* if (direction < 0) */ { 4952 if (layout_type::fast_j ()) { 4953 if (it == index2_data_.begin () + array_size_type (zero_based (*itv))) 4954 return iterator2 (*this, rank, i, j, itv, it); 4955 j = zero_based (*(it - 1)); 4956 } else { 4957 if (j == 0) 4958 return iterator2 (*this, rank, i, j, itv, it); 4959 -- j; 4960 } 4961 } 4962 } 4963 } 4964 4965 4966 class const_iterator1: 4967 public container_const_reference<coordinate_matrix>, 4968 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag, 4969 const_iterator1, value_type> { 4970 public: 4971 typedef typename coordinate_matrix::value_type value_type; 4972 typedef typename coordinate_matrix::difference_type difference_type; 4973 typedef typename coordinate_matrix::const_reference reference; 4974 typedef const typename coordinate_matrix::pointer pointer; 4975 4976 typedef const_iterator2 dual_iterator_type; 4977 typedef const_reverse_iterator2 dual_reverse_iterator_type; 4978 4979 // Construction and destruction 4980 BOOST_UBLAS_INLINE const_iterator1()4981 const_iterator1 (): 4982 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {} 4983 BOOST_UBLAS_INLINE const_iterator1(const self_type & m,int rank,size_type i,size_type j,const vector_const_subiterator_type & itv,const const_subiterator_type & it)4984 const_iterator1 (const self_type &m, int rank, size_type i, size_type j, const vector_const_subiterator_type &itv, const const_subiterator_type &it): 4985 container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {} 4986 BOOST_UBLAS_INLINE const_iterator1(const iterator1 & it)4987 const_iterator1 (const iterator1 &it): 4988 container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {} 4989 4990 // Arithmetic 4991 BOOST_UBLAS_INLINE operator ++()4992 const_iterator1 &operator ++ () { 4993 if (rank_ == 1 && layout_type::fast_i ()) 4994 ++ it_; 4995 else { 4996 i_ = index1 () + 1; 4997 if (rank_ == 1) 4998 *this = (*this) ().find1 (rank_, i_, j_, 1); 4999 } 5000 return *this; 5001 } 5002 BOOST_UBLAS_INLINE operator --()5003 const_iterator1 &operator -- () { 5004 if (rank_ == 1 && layout_type::fast_i ()) 5005 -- it_; 5006 else { 5007 i_ = index1 () - 1; 5008 if (rank_ == 1) 5009 *this = (*this) ().find1 (rank_, i_, j_, -1); 5010 } 5011 return *this; 5012 } 5013 5014 // Dereference 5015 BOOST_UBLAS_INLINE operator *() const5016 const_reference operator * () const { 5017 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ()); 5018 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ()); 5019 if (rank_ == 1) { 5020 return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()]; 5021 } else { 5022 return (*this) () (i_, j_); 5023 } 5024 } 5025 5026 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 5027 BOOST_UBLAS_INLINE 5028 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5029 typename self_type:: 5030 #endif begin() const5031 const_iterator2 begin () const { 5032 const self_type &m = (*this) (); 5033 return m.find2 (1, index1 (), 0); 5034 } 5035 BOOST_UBLAS_INLINE 5036 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5037 typename self_type:: 5038 #endif cbegin() const5039 const_iterator2 cbegin () const { 5040 return begin (); 5041 } 5042 BOOST_UBLAS_INLINE 5043 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5044 typename self_type:: 5045 #endif end() const5046 const_iterator2 end () const { 5047 const self_type &m = (*this) (); 5048 return m.find2 (1, index1 (), m.size2 ()); 5049 } 5050 BOOST_UBLAS_INLINE 5051 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5052 typename self_type:: 5053 #endif cend() const5054 const_iterator2 cend () const { 5055 return end (); 5056 } 5057 BOOST_UBLAS_INLINE 5058 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5059 typename self_type:: 5060 #endif rbegin() const5061 const_reverse_iterator2 rbegin () const { 5062 return const_reverse_iterator2 (end ()); 5063 } 5064 BOOST_UBLAS_INLINE 5065 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5066 typename self_type:: 5067 #endif crbegin() const5068 const_reverse_iterator2 crbegin () const { 5069 return rbegin (); 5070 } 5071 BOOST_UBLAS_INLINE 5072 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5073 typename self_type:: 5074 #endif rend() const5075 const_reverse_iterator2 rend () const { 5076 return const_reverse_iterator2 (begin ()); 5077 } 5078 BOOST_UBLAS_INLINE 5079 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5080 typename self_type:: 5081 #endif crend() const5082 const_reverse_iterator2 crend () const { 5083 return rend (); 5084 } 5085 #endif 5086 5087 // Indices 5088 BOOST_UBLAS_INLINE index1() const5089 size_type index1 () const { 5090 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ()); 5091 if (rank_ == 1) { 5092 BOOST_UBLAS_CHECK (layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ()); 5093 return layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)); 5094 } else { 5095 return i_; 5096 } 5097 } 5098 BOOST_UBLAS_INLINE index2() const5099 size_type index2 () const { 5100 if (rank_ == 1) { 5101 BOOST_UBLAS_CHECK (layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ()); 5102 return layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)); 5103 } else { 5104 return j_; 5105 } 5106 } 5107 5108 // Assignment 5109 BOOST_UBLAS_INLINE operator =(const const_iterator1 & it)5110 const_iterator1 &operator = (const const_iterator1 &it) { 5111 container_const_reference<self_type>::assign (&it ()); 5112 rank_ = it.rank_; 5113 i_ = it.i_; 5114 j_ = it.j_; 5115 itv_ = it.itv_; 5116 it_ = it.it_; 5117 return *this; 5118 } 5119 5120 // Comparison 5121 BOOST_UBLAS_INLINE operator ==(const const_iterator1 & it) const5122 bool operator == (const const_iterator1 &it) const { 5123 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 5124 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ()); 5125 if (rank_ == 1 || it.rank_ == 1) { 5126 return it_ == it.it_; 5127 } else { 5128 return i_ == it.i_ && j_ == it.j_; 5129 } 5130 } 5131 5132 private: 5133 int rank_; 5134 size_type i_; 5135 size_type j_; 5136 vector_const_subiterator_type itv_; 5137 const_subiterator_type it_; 5138 }; 5139 5140 BOOST_UBLAS_INLINE begin1() const5141 const_iterator1 begin1 () const { 5142 return find1 (0, 0, 0); 5143 } 5144 BOOST_UBLAS_INLINE cbegin1() const5145 const_iterator1 cbegin1 () const { 5146 return begin1 (); 5147 } 5148 BOOST_UBLAS_INLINE end1() const5149 const_iterator1 end1 () const { 5150 return find1 (0, size1_, 0); 5151 } 5152 BOOST_UBLAS_INLINE cend1() const5153 const_iterator1 cend1 () const { 5154 return end1 (); 5155 } 5156 5157 class iterator1: 5158 public container_reference<coordinate_matrix>, 5159 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag, 5160 iterator1, value_type> { 5161 public: 5162 typedef typename coordinate_matrix::value_type value_type; 5163 typedef typename coordinate_matrix::difference_type difference_type; 5164 typedef typename coordinate_matrix::true_reference reference; 5165 typedef typename coordinate_matrix::pointer pointer; 5166 5167 typedef iterator2 dual_iterator_type; 5168 typedef reverse_iterator2 dual_reverse_iterator_type; 5169 5170 // Construction and destruction 5171 BOOST_UBLAS_INLINE iterator1()5172 iterator1 (): 5173 container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {} 5174 BOOST_UBLAS_INLINE iterator1(self_type & m,int rank,size_type i,size_type j,const vector_subiterator_type & itv,const subiterator_type & it)5175 iterator1 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it): 5176 container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {} 5177 5178 // Arithmetic 5179 BOOST_UBLAS_INLINE operator ++()5180 iterator1 &operator ++ () { 5181 if (rank_ == 1 && layout_type::fast_i ()) 5182 ++ it_; 5183 else { 5184 i_ = index1 () + 1; 5185 if (rank_ == 1) 5186 *this = (*this) ().find1 (rank_, i_, j_, 1); 5187 } 5188 return *this; 5189 } 5190 BOOST_UBLAS_INLINE operator --()5191 iterator1 &operator -- () { 5192 if (rank_ == 1 && layout_type::fast_i ()) 5193 -- it_; 5194 else { 5195 i_ = index1 () - 1; 5196 if (rank_ == 1) 5197 *this = (*this) ().find1 (rank_, i_, j_, -1); 5198 } 5199 return *this; 5200 } 5201 5202 // Dereference 5203 BOOST_UBLAS_INLINE operator *() const5204 reference operator * () const { 5205 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ()); 5206 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ()); 5207 if (rank_ == 1) { 5208 return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()]; 5209 } else { 5210 return (*this) ().at_element (i_, j_); 5211 } 5212 } 5213 5214 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 5215 BOOST_UBLAS_INLINE 5216 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5217 typename self_type:: 5218 #endif begin() const5219 iterator2 begin () const { 5220 self_type &m = (*this) (); 5221 return m.find2 (1, index1 (), 0); 5222 } 5223 BOOST_UBLAS_INLINE 5224 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5225 typename self_type:: 5226 #endif end() const5227 iterator2 end () const { 5228 self_type &m = (*this) (); 5229 return m.find2 (1, index1 (), m.size2 ()); 5230 } 5231 BOOST_UBLAS_INLINE 5232 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5233 typename self_type:: 5234 #endif rbegin() const5235 reverse_iterator2 rbegin () const { 5236 return reverse_iterator2 (end ()); 5237 } 5238 BOOST_UBLAS_INLINE 5239 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5240 typename self_type:: 5241 #endif rend() const5242 reverse_iterator2 rend () const { 5243 return reverse_iterator2 (begin ()); 5244 } 5245 #endif 5246 5247 // Indices 5248 BOOST_UBLAS_INLINE index1() const5249 size_type index1 () const { 5250 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ()); 5251 if (rank_ == 1) { 5252 BOOST_UBLAS_CHECK (layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ()); 5253 return layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)); 5254 } else { 5255 return i_; 5256 } 5257 } 5258 BOOST_UBLAS_INLINE index2() const5259 size_type index2 () const { 5260 if (rank_ == 1) { 5261 BOOST_UBLAS_CHECK (layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ()); 5262 return layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)); 5263 } else { 5264 return j_; 5265 } 5266 } 5267 5268 // Assignment 5269 BOOST_UBLAS_INLINE operator =(const iterator1 & it)5270 iterator1 &operator = (const iterator1 &it) { 5271 container_reference<self_type>::assign (&it ()); 5272 rank_ = it.rank_; 5273 i_ = it.i_; 5274 j_ = it.j_; 5275 itv_ = it.itv_; 5276 it_ = it.it_; 5277 return *this; 5278 } 5279 5280 // Comparison 5281 BOOST_UBLAS_INLINE operator ==(const iterator1 & it) const5282 bool operator == (const iterator1 &it) const { 5283 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 5284 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ()); 5285 if (rank_ == 1 || it.rank_ == 1) { 5286 return it_ == it.it_; 5287 } else { 5288 return i_ == it.i_ && j_ == it.j_; 5289 } 5290 } 5291 5292 private: 5293 int rank_; 5294 size_type i_; 5295 size_type j_; 5296 vector_subiterator_type itv_; 5297 subiterator_type it_; 5298 5299 friend class const_iterator1; 5300 }; 5301 5302 BOOST_UBLAS_INLINE begin1()5303 iterator1 begin1 () { 5304 return find1 (0, 0, 0); 5305 } 5306 BOOST_UBLAS_INLINE end1()5307 iterator1 end1 () { 5308 return find1 (0, size1_, 0); 5309 } 5310 5311 class const_iterator2: 5312 public container_const_reference<coordinate_matrix>, 5313 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag, 5314 const_iterator2, value_type> { 5315 public: 5316 typedef typename coordinate_matrix::value_type value_type; 5317 typedef typename coordinate_matrix::difference_type difference_type; 5318 typedef typename coordinate_matrix::const_reference reference; 5319 typedef const typename coordinate_matrix::pointer pointer; 5320 5321 typedef const_iterator1 dual_iterator_type; 5322 typedef const_reverse_iterator1 dual_reverse_iterator_type; 5323 5324 // Construction and destruction 5325 BOOST_UBLAS_INLINE const_iterator2()5326 const_iterator2 (): 5327 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {} 5328 BOOST_UBLAS_INLINE const_iterator2(const self_type & m,int rank,size_type i,size_type j,const vector_const_subiterator_type itv,const const_subiterator_type & it)5329 const_iterator2 (const self_type &m, int rank, size_type i, size_type j, const vector_const_subiterator_type itv, const const_subiterator_type &it): 5330 container_const_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {} 5331 BOOST_UBLAS_INLINE const_iterator2(const iterator2 & it)5332 const_iterator2 (const iterator2 &it): 5333 container_const_reference<self_type> (it ()), rank_ (it.rank_), i_ (it.i_), j_ (it.j_), itv_ (it.itv_), it_ (it.it_) {} 5334 5335 // Arithmetic 5336 BOOST_UBLAS_INLINE operator ++()5337 const_iterator2 &operator ++ () { 5338 if (rank_ == 1 && layout_type::fast_j ()) 5339 ++ it_; 5340 else { 5341 j_ = index2 () + 1; 5342 if (rank_ == 1) 5343 *this = (*this) ().find2 (rank_, i_, j_, 1); 5344 } 5345 return *this; 5346 } 5347 BOOST_UBLAS_INLINE operator --()5348 const_iterator2 &operator -- () { 5349 if (rank_ == 1 && layout_type::fast_j ()) 5350 -- it_; 5351 else { 5352 j_ = index2 () - 1; 5353 if (rank_ == 1) 5354 *this = (*this) ().find2 (rank_, i_, j_, -1); 5355 } 5356 return *this; 5357 } 5358 5359 // Dereference 5360 BOOST_UBLAS_INLINE operator *() const5361 const_reference operator * () const { 5362 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ()); 5363 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ()); 5364 if (rank_ == 1) { 5365 return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()]; 5366 } else { 5367 return (*this) () (i_, j_); 5368 } 5369 } 5370 5371 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 5372 BOOST_UBLAS_INLINE 5373 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5374 typename self_type:: 5375 #endif begin() const5376 const_iterator1 begin () const { 5377 const self_type &m = (*this) (); 5378 return m.find1 (1, 0, index2 ()); 5379 } 5380 BOOST_UBLAS_INLINE 5381 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5382 typename self_type:: 5383 #endif cbegin() const5384 const_iterator1 cbegin () const { 5385 return begin (); 5386 } 5387 BOOST_UBLAS_INLINE 5388 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5389 typename self_type:: 5390 #endif end() const5391 const_iterator1 end () const { 5392 const self_type &m = (*this) (); 5393 return m.find1 (1, m.size1 (), index2 ()); 5394 } 5395 BOOST_UBLAS_INLINE 5396 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5397 typename self_type:: 5398 #endif cend() const5399 const_iterator1 cend () const { 5400 return end (); 5401 } 5402 BOOST_UBLAS_INLINE 5403 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5404 typename self_type:: 5405 #endif rbegin() const5406 const_reverse_iterator1 rbegin () const { 5407 return const_reverse_iterator1 (end ()); 5408 } 5409 BOOST_UBLAS_INLINE 5410 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5411 typename self_type:: 5412 #endif crbegin() const5413 const_reverse_iterator1 crbegin () const { 5414 return rbegin (); 5415 } 5416 BOOST_UBLAS_INLINE 5417 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5418 typename self_type:: 5419 #endif rend() const5420 const_reverse_iterator1 rend () const { 5421 return const_reverse_iterator1 (begin ()); 5422 } 5423 BOOST_UBLAS_INLINE 5424 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5425 typename self_type:: 5426 #endif crend() const5427 const_reverse_iterator1 crend () const { 5428 return rend (); 5429 } 5430 #endif 5431 5432 // Indices 5433 BOOST_UBLAS_INLINE index1() const5434 size_type index1 () const { 5435 if (rank_ == 1) { 5436 BOOST_UBLAS_CHECK (layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ()); 5437 return layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)); 5438 } else { 5439 return i_; 5440 } 5441 } 5442 BOOST_UBLAS_INLINE index2() const5443 size_type index2 () const { 5444 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ()); 5445 if (rank_ == 1) { 5446 BOOST_UBLAS_CHECK (layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ()); 5447 return layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)); 5448 } else { 5449 return j_; 5450 } 5451 } 5452 5453 // Assignment 5454 BOOST_UBLAS_INLINE operator =(const const_iterator2 & it)5455 const_iterator2 &operator = (const const_iterator2 &it) { 5456 container_const_reference<self_type>::assign (&it ()); 5457 rank_ = it.rank_; 5458 i_ = it.i_; 5459 j_ = it.j_; 5460 itv_ = it.itv_; 5461 it_ = it.it_; 5462 return *this; 5463 } 5464 5465 // Comparison 5466 BOOST_UBLAS_INLINE operator ==(const const_iterator2 & it) const5467 bool operator == (const const_iterator2 &it) const { 5468 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 5469 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ()); 5470 if (rank_ == 1 || it.rank_ == 1) { 5471 return it_ == it.it_; 5472 } else { 5473 return i_ == it.i_ && j_ == it.j_; 5474 } 5475 } 5476 5477 private: 5478 int rank_; 5479 size_type i_; 5480 size_type j_; 5481 vector_const_subiterator_type itv_; 5482 const_subiterator_type it_; 5483 }; 5484 5485 BOOST_UBLAS_INLINE begin2() const5486 const_iterator2 begin2 () const { 5487 return find2 (0, 0, 0); 5488 } 5489 BOOST_UBLAS_INLINE cbegin2() const5490 const_iterator2 cbegin2 () const { 5491 return begin2 (); 5492 } 5493 BOOST_UBLAS_INLINE end2() const5494 const_iterator2 end2 () const { 5495 return find2 (0, 0, size2_); 5496 } 5497 BOOST_UBLAS_INLINE cend2() const5498 const_iterator2 cend2 () const { 5499 return end2 (); 5500 } 5501 5502 class iterator2: 5503 public container_reference<coordinate_matrix>, 5504 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag, 5505 iterator2, value_type> { 5506 public: 5507 typedef typename coordinate_matrix::value_type value_type; 5508 typedef typename coordinate_matrix::difference_type difference_type; 5509 typedef typename coordinate_matrix::true_reference reference; 5510 typedef typename coordinate_matrix::pointer pointer; 5511 5512 typedef iterator1 dual_iterator_type; 5513 typedef reverse_iterator1 dual_reverse_iterator_type; 5514 5515 // Construction and destruction 5516 BOOST_UBLAS_INLINE iterator2()5517 iterator2 (): 5518 container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), it_ () {} 5519 BOOST_UBLAS_INLINE iterator2(self_type & m,int rank,size_type i,size_type j,const vector_subiterator_type & itv,const subiterator_type & it)5520 iterator2 (self_type &m, int rank, size_type i, size_type j, const vector_subiterator_type &itv, const subiterator_type &it): 5521 container_reference<self_type> (m), rank_ (rank), i_ (i), j_ (j), itv_ (itv), it_ (it) {} 5522 5523 // Arithmetic 5524 BOOST_UBLAS_INLINE operator ++()5525 iterator2 &operator ++ () { 5526 if (rank_ == 1 && layout_type::fast_j ()) 5527 ++ it_; 5528 else { 5529 j_ = index2 () + 1; 5530 if (rank_ == 1) 5531 *this = (*this) ().find2 (rank_, i_, j_, 1); 5532 } 5533 return *this; 5534 } 5535 BOOST_UBLAS_INLINE operator --()5536 iterator2 &operator -- () { 5537 if (rank_ == 1 && layout_type::fast_j ()) 5538 -- it_; 5539 else { 5540 j_ = index2 (); 5541 if (rank_ == 1) 5542 *this = (*this) ().find2 (rank_, i_, j_, -1); 5543 } 5544 return *this; 5545 } 5546 5547 // Dereference 5548 BOOST_UBLAS_INLINE operator *() const5549 reference operator * () const { 5550 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ()); 5551 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ()); 5552 if (rank_ == 1) { 5553 return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()]; 5554 } else { 5555 return (*this) ().at_element (i_, j_); 5556 } 5557 } 5558 5559 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION 5560 BOOST_UBLAS_INLINE 5561 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5562 typename self_type:: 5563 #endif begin() const5564 iterator1 begin () const { 5565 self_type &m = (*this) (); 5566 return m.find1 (1, 0, index2 ()); 5567 } 5568 BOOST_UBLAS_INLINE 5569 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5570 typename self_type:: 5571 #endif end() const5572 iterator1 end () const { 5573 self_type &m = (*this) (); 5574 return m.find1 (1, m.size1 (), index2 ()); 5575 } 5576 BOOST_UBLAS_INLINE 5577 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5578 typename self_type:: 5579 #endif rbegin() const5580 reverse_iterator1 rbegin () const { 5581 return reverse_iterator1 (end ()); 5582 } 5583 BOOST_UBLAS_INLINE 5584 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION 5585 typename self_type:: 5586 #endif rend() const5587 reverse_iterator1 rend () const { 5588 return reverse_iterator1 (begin ()); 5589 } 5590 #endif 5591 5592 // Indices 5593 BOOST_UBLAS_INLINE index1() const5594 size_type index1 () const { 5595 if (rank_ == 1) { 5596 BOOST_UBLAS_CHECK (layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size1 (), bad_index ()); 5597 return layout_type::index_M ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)); 5598 } else { 5599 return i_; 5600 } 5601 } 5602 BOOST_UBLAS_INLINE index2() const5603 size_type index2 () const { 5604 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ()); 5605 if (rank_ == 1) { 5606 BOOST_UBLAS_CHECK (layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)) < (*this) ().size2 (), bad_index ()); 5607 return layout_type::index_m ((*this) ().zero_based (*itv_), (*this) ().zero_based (*it_)); 5608 } else { 5609 return j_; 5610 } 5611 } 5612 5613 // Assignment 5614 BOOST_UBLAS_INLINE operator =(const iterator2 & it)5615 iterator2 &operator = (const iterator2 &it) { 5616 container_reference<self_type>::assign (&it ()); 5617 rank_ = it.rank_; 5618 i_ = it.i_; 5619 j_ = it.j_; 5620 itv_ = it.itv_; 5621 it_ = it.it_; 5622 return *this; 5623 } 5624 5625 // Comparison 5626 BOOST_UBLAS_INLINE operator ==(const iterator2 & it) const5627 bool operator == (const iterator2 &it) const { 5628 BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ()); 5629 // BOOST_UBLAS_CHECK (rank_ == it.rank_, internal_logic ()); 5630 if (rank_ == 1 || it.rank_ == 1) { 5631 return it_ == it.it_; 5632 } else { 5633 return i_ == it.i_ && j_ == it.j_; 5634 } 5635 } 5636 5637 private: 5638 int rank_; 5639 size_type i_; 5640 size_type j_; 5641 vector_subiterator_type itv_; 5642 subiterator_type it_; 5643 5644 friend class const_iterator2; 5645 }; 5646 5647 BOOST_UBLAS_INLINE begin2()5648 iterator2 begin2 () { 5649 return find2 (0, 0, 0); 5650 } 5651 BOOST_UBLAS_INLINE end2()5652 iterator2 end2 () { 5653 return find2 (0, 0, size2_); 5654 } 5655 5656 // Reverse iterators 5657 5658 BOOST_UBLAS_INLINE rbegin1() const5659 const_reverse_iterator1 rbegin1 () const { 5660 return const_reverse_iterator1 (end1 ()); 5661 } 5662 BOOST_UBLAS_INLINE crbegin1() const5663 const_reverse_iterator1 crbegin1 () const { 5664 return rbegin1 (); 5665 } 5666 BOOST_UBLAS_INLINE rend1() const5667 const_reverse_iterator1 rend1 () const { 5668 return const_reverse_iterator1 (begin1 ()); 5669 } 5670 BOOST_UBLAS_INLINE crend1() const5671 const_reverse_iterator1 crend1 () const { 5672 return rend1 (); 5673 } 5674 5675 BOOST_UBLAS_INLINE rbegin1()5676 reverse_iterator1 rbegin1 () { 5677 return reverse_iterator1 (end1 ()); 5678 } 5679 BOOST_UBLAS_INLINE rend1()5680 reverse_iterator1 rend1 () { 5681 return reverse_iterator1 (begin1 ()); 5682 } 5683 5684 BOOST_UBLAS_INLINE rbegin2() const5685 const_reverse_iterator2 rbegin2 () const { 5686 return const_reverse_iterator2 (end2 ()); 5687 } 5688 BOOST_UBLAS_INLINE crbegin2() const5689 const_reverse_iterator2 crbegin2 () const { 5690 return rbegin2 (); 5691 } 5692 BOOST_UBLAS_INLINE rend2() const5693 const_reverse_iterator2 rend2 () const { 5694 return const_reverse_iterator2 (begin2 ()); 5695 } 5696 BOOST_UBLAS_INLINE crend2() const5697 const_reverse_iterator2 crend2 () const { 5698 return rend2 (); 5699 } 5700 5701 BOOST_UBLAS_INLINE rbegin2()5702 reverse_iterator2 rbegin2 () { 5703 return reverse_iterator2 (end2 ()); 5704 } 5705 BOOST_UBLAS_INLINE rend2()5706 reverse_iterator2 rend2 () { 5707 return reverse_iterator2 (begin2 ()); 5708 } 5709 5710 // Serialization 5711 template<class Archive> serialize(Archive & ar,const unsigned int)5712 void serialize(Archive & ar, const unsigned int /* file_version */){ 5713 serialization::collection_size_type s1 (size1_); 5714 serialization::collection_size_type s2 (size2_); 5715 ar & serialization::make_nvp("size1",s1); 5716 ar & serialization::make_nvp("size2",s2); 5717 if (Archive::is_loading::value) { 5718 size1_ = s1; 5719 size2_ = s2; 5720 } 5721 ar & serialization::make_nvp("capacity", capacity_); 5722 ar & serialization::make_nvp("filled", filled_); 5723 ar & serialization::make_nvp("sorted_filled", sorted_filled_); 5724 ar & serialization::make_nvp("sorted", sorted_); 5725 ar & serialization::make_nvp("index1_data", index1_data_); 5726 ar & serialization::make_nvp("index2_data", index2_data_); 5727 ar & serialization::make_nvp("value_data", value_data_); 5728 storage_invariants(); 5729 } 5730 5731 private: storage_invariants() const5732 void storage_invariants () const 5733 { 5734 BOOST_UBLAS_CHECK (capacity_ == index1_data_.size (), internal_logic ()); 5735 BOOST_UBLAS_CHECK (capacity_ == index2_data_.size (), internal_logic ()); 5736 BOOST_UBLAS_CHECK (capacity_ == value_data_.size (), internal_logic ()); 5737 BOOST_UBLAS_CHECK (filled_ <= capacity_, internal_logic ()); 5738 BOOST_UBLAS_CHECK (sorted_filled_ <= filled_, internal_logic ()); 5739 BOOST_UBLAS_CHECK (sorted_ == (sorted_filled_ == filled_), internal_logic ()); 5740 } 5741 5742 size_type size1_; 5743 size_type size2_; 5744 array_size_type capacity_; 5745 mutable array_size_type filled_; 5746 mutable array_size_type sorted_filled_; 5747 mutable bool sorted_; 5748 mutable index_array_type index1_data_; 5749 mutable index_array_type index2_data_; 5750 mutable value_array_type value_data_; 5751 static const value_type zero_; 5752 5753 BOOST_UBLAS_INLINE zero_based(size_type k_based_index)5754 static size_type zero_based (size_type k_based_index) { 5755 return k_based_index - IB; 5756 } 5757 BOOST_UBLAS_INLINE k_based(size_type zero_based_index)5758 static size_type k_based (size_type zero_based_index) { 5759 return zero_based_index + IB; 5760 } 5761 5762 friend class iterator1; 5763 friend class iterator2; 5764 friend class const_iterator1; 5765 friend class const_iterator2; 5766 }; 5767 5768 template<class T, class L, std::size_t IB, class IA, class TA> 5769 const typename coordinate_matrix<T, L, IB, IA, TA>::value_type coordinate_matrix<T, L, IB, IA, TA>::zero_ = value_type/*zero*/(); 5770 5771 }}} 5772 5773 #endif 5774