2 // Copyright (c) 2000-2007
3 // Joerg Walter, Mathias Koch, Gunter Winkler
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)
9 // The authors gratefully acknowledge the support of
10 // GeNeSys mbH & Co. KG in producing this work.
13 #ifndef _BOOST_UBLAS_MATRIX_SPARSE_
14 #define _BOOST_UBLAS_MATRIX_SPARSE_
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>
23 // Iterators based on ideas of Jeremy Siek
25 namespace boost { namespace numeric { namespace ublas {
27 #ifdef BOOST_UBLAS_STRICT_MATRIX_SPARSE
30 class sparse_matrix_element:
31 public container_reference<M> {
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;
41 // Proxied element operations
43 const_pointer p = (*this) ().find_element (i_, j_);
47 d_ = value_type/*zero*/();
50 void set (const value_type &s) const {
51 pointer p = (*this) ().find_element (i_, j_);
53 (*this) ().insert_element (i_, j_, s);
59 // Construction and destruction
61 sparse_matrix_element (matrix_type &m, size_type i, size_type j):
62 container_reference<matrix_type> (m), i_ (i), j_ (j) {
65 sparse_matrix_element (const sparse_matrix_element &p):
66 container_reference<matrix_type> (p), i_ (p.i_), j_ (p.j_) {}
68 ~sparse_matrix_element () {
73 sparse_matrix_element &operator = (const sparse_matrix_element &p) {
74 // Overide the implict copy assignment
81 sparse_matrix_element &operator = (const D &d) {
87 sparse_matrix_element &operator += (const D &d) {
95 sparse_matrix_element &operator -= (const D &d) {
103 sparse_matrix_element &operator *= (const D &d) {
111 sparse_matrix_element &operator /= (const D &d) {
121 bool operator == (const D &d) const {
127 bool operator != (const D &d) const {
132 // Conversion - weak link in proxy as d_ is not a perfect alias for the element
134 operator const_reference () const {
139 // Conversion to reference - may be invalidated
141 value_type& ref () const {
142 const pointer p = (*this) ().find_element (i_, j_);
144 return (*this) ().insert_element (i_, j_, value_type/*zero*/());
152 mutable value_type d_;
156 * Generalise explicit reference access
160 struct element_reference<sparse_matrix_element<V> > {
161 typedef typename V::value_type& reference;
162 static reference get_reference (const sparse_matrix_element<V>& sve)
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;
180 static const unsigned plus_complexity = type_traits<element_type>::plus_complexity;
181 static const unsigned multiplies_complexity = type_traits<element_type>::multiplies_complexity;
185 real_type real (const_reference t) {
186 return type_traits<element_type>::real (t);
190 real_type imag (const_reference t) {
191 return type_traits<element_type>::imag (t);
195 value_type conj (const_reference t) {
196 return type_traits<element_type>::conj (t);
201 real_type type_abs (const_reference t) {
202 return type_traits<element_type>::type_abs (t);
206 value_type type_sqrt (const_reference t) {
207 return type_traits<element_type>::type_sqrt (t);
212 real_type norm_1 (const_reference t) {
213 return type_traits<element_type>::norm_1 (t);
217 real_type norm_2 (const_reference t) {
218 return type_traits<element_type>::norm_2 (t);
222 real_type norm_inf (const_reference t) {
223 return type_traits<element_type>::norm_inf (t);
228 bool equals (const_reference t1, const_reference t2) {
229 return type_traits<element_type>::equals (t1, t2);
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;
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;
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;
249 /** \brief Index map based sparse matrix of values of type \c T
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.
262 * \sa fwd.hpp, storage_sparse.hpp
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
267 template<class T, class L, class A>
269 public matrix_container<mapped_matrix<T, L, A> > {
271 typedef T &true_reference;
273 typedef const T * const_pointer;
274 typedef L layout_type;
275 typedef mapped_matrix<T, L, A> self_type;
277 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
278 using matrix_container<self_type>::operator ();
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;
288 typedef sparse_matrix_element<self_type> reference;
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;
297 // Construction and destruction
300 matrix_container<self_type> (),
301 size1_ (0), size2_ (0), data_ () {}
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));
309 mapped_matrix (const mapped_matrix &m):
310 matrix_container<self_type> (),
311 size1_ (m.size1_), size2_ (m.size2_), data_ (m.data_) {}
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);
323 size_type size1 () const {
327 size_type size2 () const {
331 size_type nnz_capacity () const {
332 return detail::map_capacity (data ());
335 size_type nnz () const {
336 return data (). size ();
341 const array_type &data () const {
345 array_type &data () {
352 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_;
361 void resize (size_type size1, size_type size2, bool preserve = true) {
362 // FIXME preserve unimplemented
363 BOOST_UBLAS_CHECK (!preserve, internal_logic ());
371 void reserve (size_type non_zeros, bool preserve = true) {
372 detail::map_reserve (data (), restrict_capacity (non_zeros));
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));
381 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 ())
386 BOOST_UBLAS_CHECK ((*it).first == element, internal_logic ()); // broken map
387 return &(*it).second;
392 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 ())
397 BOOST_UBLAS_CHECK ((*it).first == element, internal_logic ()); // broken map
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;
408 return reference (*this, i, j);
412 // Element assingment
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;
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 ())
439 mapped_matrix &operator = (const mapped_matrix &m) {
447 template<class C> // Container assignment without temporary
449 mapped_matrix &operator = (const matrix_container<C> &m) {
450 resize (m ().size1 (), m ().size2 (), false);
455 mapped_matrix &assign_temporary (mapped_matrix &m) {
461 mapped_matrix &operator = (const matrix_expression<AE> &ae) {
462 self_type temporary (ae, detail::map_capacity (data ()));
463 return assign_temporary (temporary);
467 mapped_matrix &assign (const matrix_expression<AE> &ae) {
468 matrix_assign<scalar_assign> (*this, 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);
477 template<class C> // Container assignment without temporary
479 mapped_matrix &operator += (const matrix_container<C> &m) {
485 mapped_matrix &plus_assign (const matrix_expression<AE> &ae) {
486 matrix_assign<scalar_plus_assign> (*this, 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);
495 template<class C> // Container assignment without temporary
497 mapped_matrix &operator -= (const matrix_container<C> &m) {
503 mapped_matrix &minus_assign (const matrix_expression<AE> &ae) {
504 matrix_assign<scalar_minus_assign> (*this, ae);
509 mapped_matrix& operator *= (const AT &at) {
510 matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
515 mapped_matrix& operator /= (const AT &at) {
516 matrix_assign_scalar<scalar_divides_assign> (*this, at);
522 void swap (mapped_matrix &m) {
524 std::swap (size1_, m.size1_);
525 std::swap (size2_, m.size2_);
526 data ().swap (m.data ());
530 friend void swap (mapped_matrix &m1, mapped_matrix &m2) {
536 // Use storage iterator
537 typedef typename A::const_iterator const_subiterator_type;
538 typedef typename A::iterator subiterator_type;
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
550 class const_iterator1;
552 class const_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;
560 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
561 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_);
570 if ((index1 >= i && index2 == j) || (i >= size1_))
573 } else /* if (direction < 0) */ {
574 if ((index1 <= i && index2 == j) || (i == 0))
578 it = data ().lower_bound (layout_type::address (i, size1_, j, size2_));
580 if (rank == 1 && index2 != j) {
583 else /* if (direction < 0) */
587 return const_iterator1 (*this, rank, i, j, it);
589 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
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_);
599 if ((index1 >= i && index2 == j) || (i >= size1_))
602 } else /* if (direction < 0) */ {
603 if ((index1 <= i && index2 == j) || (i == 0))
607 it = data ().lower_bound (layout_type::address (i, size1_, j, size2_));
609 if (rank == 1 && index2 != j) {
612 else /* if (direction < 0) */
616 return iterator1 (*this, rank, i, j, it);
618 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
619 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_);
628 if ((index2 >= j && index1 == i) || (j >= size2_))
631 } else /* if (direction < 0) */ {
632 if ((index2 <= j && index1 == i) || (j == 0))
636 it = data ().lower_bound (layout_type::address (i, size1_, j, size2_));
638 if (rank == 1 && index1 != i) {
641 else /* if (direction < 0) */
645 return const_iterator2 (*this, rank, i, j, it);
647 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
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_);
657 if ((index2 >= j && index1 == i) || (j >= size2_))
660 } else /* if (direction < 0) */ {
661 if ((index2 <= j && index1 == i) || (j == 0))
665 it = data ().lower_bound (layout_type::address (i, size1_, j, size2_));
667 if (rank == 1 && index1 != i) {
670 else /* if (direction < 0) */
674 return iterator2 (*this, rank, i, j, it);
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> {
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;
688 typedef const_iterator2 dual_iterator_type;
689 typedef const_reverse_iterator2 dual_reverse_iterator_type;
691 // Construction and destruction
694 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), 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) {}
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_) {}
704 const_iterator1 &operator ++ () {
705 if (rank_ == 1 && layout_type::fast_i ())
708 *this = (*this) ().find1 (rank_, index1 () + 1, j_, 1);
712 const_iterator1 &operator -- () {
713 if (rank_ == 1 && layout_type::fast_i ())
716 *this = (*this) ().find1 (rank_, index1 () - 1, j_, -1);
722 const_reference operator * () const {
723 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
724 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
726 return (*it_).second;
728 return (*this) () (i_, j_);
732 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
734 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
737 const_iterator2 begin () const {
738 const self_type &m = (*this) ();
739 return m.find2 (1, index1 (), 0);
742 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
745 const_iterator2 cbegin () const {
749 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
752 const_iterator2 end () const {
753 const self_type &m = (*this) ();
754 return m.find2 (1, index1 (), m.size2 ());
757 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
760 const_iterator2 cend () const {
764 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
767 const_reverse_iterator2 rbegin () const {
768 return const_reverse_iterator2 (end ());
771 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
774 const_reverse_iterator2 crbegin () const {
778 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
781 const_reverse_iterator2 rend () const {
782 return const_reverse_iterator2 (begin ());
785 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
788 const_reverse_iterator2 crend () const {
795 size_type index1 () const {
796 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
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 ());
806 size_type index2 () const {
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 ());
818 const_iterator1 &operator = (const const_iterator1 &it) {
819 container_const_reference<self_type>::assign (&it ());
829 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_;
835 return i_ == it.i_ && j_ == it.j_;
843 const_subiterator_type it_;
847 const_iterator1 begin1 () const {
848 return find1 (0, 0, 0);
851 const_iterator1 cbegin1 () const {
855 const_iterator1 end1 () const {
856 return find1 (0, size1_, 0);
859 const_iterator1 cend1 () const {
864 public container_reference<mapped_matrix>,
865 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
866 iterator1, value_type> {
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;
873 typedef iterator2 dual_iterator_type;
874 typedef reverse_iterator2 dual_reverse_iterator_type;
876 // Construction and destruction
879 container_reference<self_type> (), rank_ (), i_ (), j_ (), 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) {}
886 iterator1 &operator ++ () {
887 if (rank_ == 1 && layout_type::fast_i ())
890 *this = (*this) ().find1 (rank_, index1 () + 1, j_, 1);
894 iterator1 &operator -- () {
895 if (rank_ == 1 && layout_type::fast_i ())
898 *this = (*this) ().find1 (rank_, index1 () - 1, j_, -1);
904 reference operator * () const {
905 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
906 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
908 return (*it_).second;
910 return (*this) ().at_element (i_, j_);
914 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
916 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
919 iterator2 begin () const {
920 self_type &m = (*this) ();
921 return m.find2 (1, index1 (), 0);
924 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
927 iterator2 end () const {
928 self_type &m = (*this) ();
929 return m.find2 (1, index1 (), m.size2 ());
932 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
935 reverse_iterator2 rbegin () const {
936 return reverse_iterator2 (end ());
939 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
942 reverse_iterator2 rend () const {
943 return reverse_iterator2 (begin ());
949 size_type index1 () const {
950 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
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 ());
960 size_type index2 () const {
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 ());
972 iterator1 &operator = (const iterator1 &it) {
973 container_reference<self_type>::assign (&it ());
983 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_;
989 return i_ == it.i_ && j_ == it.j_;
997 subiterator_type it_;
999 friend class const_iterator1;
1003 iterator1 begin1 () {
1004 return find1 (0, 0, 0);
1008 return find1 (0, size1_, 0);
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> {
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;
1021 typedef const_iterator1 dual_iterator_type;
1022 typedef const_reverse_iterator1 dual_reverse_iterator_type;
1024 // Construction and destruction
1027 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), 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) {}
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_) {}
1037 const_iterator2 &operator ++ () {
1038 if (rank_ == 1 && layout_type::fast_j ())
1041 *this = (*this) ().find2 (rank_, i_, index2 () + 1, 1);
1045 const_iterator2 &operator -- () {
1046 if (rank_ == 1 && layout_type::fast_j ())
1049 *this = (*this) ().find2 (rank_, i_, index2 () - 1, -1);
1055 const_reference operator * () const {
1056 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
1057 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
1059 return (*it_).second;
1061 return (*this) () (i_, j_);
1065 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1067 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1068 typename self_type::
1070 const_iterator1 begin () const {
1071 const self_type &m = (*this) ();
1072 return m.find1 (1, 0, index2 ());
1075 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1076 typename self_type::
1078 const_iterator1 cbegin () const {
1082 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1083 typename self_type::
1085 const_iterator1 end () const {
1086 const self_type &m = (*this) ();
1087 return m.find1 (1, m.size1 (), index2 ());
1090 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1091 typename self_type::
1093 const_iterator1 cend () const {
1097 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1098 typename self_type::
1100 const_reverse_iterator1 rbegin () const {
1101 return const_reverse_iterator1 (end ());
1104 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1105 typename self_type::
1107 const_reverse_iterator1 crbegin () const {
1111 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1112 typename self_type::
1114 const_reverse_iterator1 rend () const {
1115 return const_reverse_iterator1 (begin ());
1118 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1119 typename self_type::
1121 const_reverse_iterator1 crend () const {
1128 size_type index1 () const {
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 ());
1138 size_type index2 () const {
1139 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
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 ());
1151 const_iterator2 &operator = (const const_iterator2 &it) {
1152 container_const_reference<self_type>::assign (&it ());
1162 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_;
1168 return i_ == it.i_ && j_ == it.j_;
1176 const_subiterator_type it_;
1180 const_iterator2 begin2 () const {
1181 return find2 (0, 0, 0);
1184 const_iterator2 cbegin2 () const {
1188 const_iterator2 end2 () const {
1189 return find2 (0, 0, size2_);
1192 const_iterator2 cend2 () const {
1197 public container_reference<mapped_matrix>,
1198 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
1199 iterator2, value_type> {
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;
1206 typedef iterator1 dual_iterator_type;
1207 typedef reverse_iterator1 dual_reverse_iterator_type;
1209 // Construction and destruction
1212 container_reference<self_type> (), rank_ (), i_ (), j_ (), 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) {}
1219 iterator2 &operator ++ () {
1220 if (rank_ == 1 && layout_type::fast_j ())
1223 *this = (*this) ().find2 (rank_, i_, index2 () + 1, 1);
1227 iterator2 &operator -- () {
1228 if (rank_ == 1 && layout_type::fast_j ())
1231 *this = (*this) ().find2 (rank_, i_, index2 () - 1, -1);
1237 reference operator * () const {
1238 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
1239 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
1241 return (*it_).second;
1243 return (*this) ().at_element (i_, j_);
1247 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
1249 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1250 typename self_type::
1252 iterator1 begin () const {
1253 self_type &m = (*this) ();
1254 return m.find1 (1, 0, index2 ());
1257 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1258 typename self_type::
1260 iterator1 end () const {
1261 self_type &m = (*this) ();
1262 return m.find1 (1, m.size1 (), index2 ());
1265 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1266 typename self_type::
1268 reverse_iterator1 rbegin () const {
1269 return reverse_iterator1 (end ());
1272 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
1273 typename self_type::
1275 reverse_iterator1 rend () const {
1276 return reverse_iterator1 (begin ());
1282 size_type index1 () const {
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 ());
1292 size_type index2 () const {
1293 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
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 ());
1305 iterator2 &operator = (const iterator2 &it) {
1306 container_reference<self_type>::assign (&it ());
1316 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_;
1322 return i_ == it.i_ && j_ == it.j_;
1330 subiterator_type it_;
1332 friend class const_iterator2;
1336 iterator2 begin2 () {
1337 return find2 (0, 0, 0);
1341 return find2 (0, 0, size2_);
1344 // Reverse iterators
1347 const_reverse_iterator1 rbegin1 () const {
1348 return const_reverse_iterator1 (end1 ());
1351 const_reverse_iterator1 crbegin1 () const {
1355 const_reverse_iterator1 rend1 () const {
1356 return const_reverse_iterator1 (begin1 ());
1359 const_reverse_iterator1 crend1 () const {
1364 reverse_iterator1 rbegin1 () {
1365 return reverse_iterator1 (end1 ());
1368 reverse_iterator1 rend1 () {
1369 return reverse_iterator1 (begin1 ());
1373 const_reverse_iterator2 rbegin2 () const {
1374 return const_reverse_iterator2 (end2 ());
1377 const_reverse_iterator2 crbegin2 () const {
1381 const_reverse_iterator2 rend2 () const {
1382 return const_reverse_iterator2 (begin2 ());
1385 const_reverse_iterator2 crend2 () const {
1390 reverse_iterator2 rbegin2 () {
1391 return reverse_iterator2 (end2 ());
1394 reverse_iterator2 rend2 () {
1395 return reverse_iterator2 (begin2 ());
1399 template<class Archive>
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) {
1409 ar & serialization::make_nvp("data", data_);
1416 static const value_type zero_;
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*/();
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> > {
1428 typedef T &true_reference;
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;
1436 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
1437 using matrix_container<self_type>::operator ();
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;
1446 typedef sparse_matrix_element<self_type> reference;
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;
1456 // Construction and destruction
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 ();
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 ();
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_) {}
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);
1484 size_type size1 () const {
1488 size_type size2 () const {
1492 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);
1499 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 ();
1506 // Storage accessors
1508 const_array_type &data () const {
1512 array_type &data () {
1518 void resize (size_type size1, size_type size2, bool preserve = true) {
1519 // FIXME preserve unimplemented
1520 BOOST_UBLAS_CHECK (!preserve, internal_logic ());
1524 data () [layout_type::size_M (size1_, size2_)] = vector_data_value_type ();
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));
1533 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 ())
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 ())
1543 BOOST_UBLAS_CHECK ((*it).first == element2, internal_logic ()); // broken map
1544 return &(*it).second;
1549 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 ())
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 ())
1559 BOOST_UBLAS_CHECK ((*itv).first == element1, internal_logic ()); // broken map
1560 return (*it).second;
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;
1572 return reference (*this, i, j);
1576 // Element assignment
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);
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;
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 ())
1595 subiterator_type it ((*itv).second.find (layout_type::index_m (i, j)));
1596 if (it == (*itv).second.end ())
1598 (*itv).second.erase (it);
1605 data_ [layout_type::size_M (size1_, size2_)] = vector_data_value_type ();
1610 mapped_vector_of_mapped_vector &operator = (const mapped_vector_of_mapped_vector &m) {
1614 data () = m.data ();
1618 template<class C> // Container assignment without temporary
1620 mapped_vector_of_mapped_vector &operator = (const matrix_container<C> &m) {
1621 resize (m ().size1 (), m ().size2 (), false);
1626 mapped_vector_of_mapped_vector &assign_temporary (mapped_vector_of_mapped_vector &m) {
1632 mapped_vector_of_mapped_vector &operator = (const matrix_expression<AE> &ae) {
1633 self_type temporary (ae);
1634 return assign_temporary (temporary);
1638 mapped_vector_of_mapped_vector &assign (const matrix_expression<AE> &ae) {
1639 matrix_assign<scalar_assign> (*this, ae);
1644 mapped_vector_of_mapped_vector& operator += (const matrix_expression<AE> &ae) {
1645 self_type temporary (*this + ae);
1646 return assign_temporary (temporary);
1648 template<class C> // Container assignment without temporary
1650 mapped_vector_of_mapped_vector &operator += (const matrix_container<C> &m) {
1656 mapped_vector_of_mapped_vector &plus_assign (const matrix_expression<AE> &ae) {
1657 matrix_assign<scalar_plus_assign> (*this, ae);
1662 mapped_vector_of_mapped_vector& operator -= (const matrix_expression<AE> &ae) {
1663 self_type temporary (*this - ae);
1664 return assign_temporary (temporary);
1666 template<class C> // Container assignment without temporary
1668 mapped_vector_of_mapped_vector &operator -= (const matrix_container<C> &m) {
1674 mapped_vector_of_mapped_vector &minus_assign (const matrix_expression<AE> &ae) {
1675 matrix_assign<scalar_minus_assign> (*this, ae);
1680 mapped_vector_of_mapped_vector& operator *= (const AT &at) {
1681 matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
1686 mapped_vector_of_mapped_vector& operator /= (const AT &at) {
1687 matrix_assign_scalar<scalar_divides_assign> (*this, at);
1693 void swap (mapped_vector_of_mapped_vector &m) {
1695 std::swap (size1_, m.size1_);
1696 std::swap (size2_, m.size2_);
1697 data ().swap (m.data ());
1701 friend void swap (mapped_vector_of_mapped_vector &m1, mapped_vector_of_mapped_vector &m2) {
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;
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
1728 class const_iterator1;
1730 class const_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;
1738 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1739 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 ());
1742 vector_const_subiterator_type itv (data ().lower_bound (layout_type::index_M (i, j)));
1743 vector_const_subiterator_type itv_end (data ().end ());
1745 return const_iterator1 (*this, rank, i, j, itv_end, (*(-- itv)).second.end ());
1747 const_subiterator_type it ((*itv).second.lower_bound (layout_type::index_m (i, j)));
1748 const_subiterator_type it_end ((*itv).second.end ());
1750 // advance to the first available major index
1751 size_type M = itv->first;
1756 m = layout_type::size_m(size1_, size2_);
1758 size_type first_i = layout_type::index_M(M,m);
1759 return const_iterator1 (*this, rank, first_i, j, itv, it);
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 ()) {
1766 return const_iterator1 (*this, rank, i, j, itv, it);
1770 return const_iterator1 (*this, rank, i, j, itv, it);
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);
1781 return const_iterator1 (*this, rank, i, j, itv, it);
1787 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1788 iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) {
1789 BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ());
1791 vector_subiterator_type itv (data ().lower_bound (layout_type::index_M (i, j)));
1792 vector_subiterator_type itv_end (data ().end ());
1794 return iterator1 (*this, rank, i, j, itv_end, (*(-- itv)).second.end ());
1796 subiterator_type it ((*itv).second.lower_bound (layout_type::index_m (i, j)));
1797 subiterator_type it_end ((*itv).second.end ());
1799 // advance to the first available major index
1800 size_type M = itv->first;
1805 m = layout_type::size_m(size1_, size2_);
1807 size_type first_i = layout_type::index_M(M,m);
1808 return iterator1 (*this, rank, first_i, j, itv, it);
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 ()) {
1815 return iterator1 (*this, rank, i, j, itv, it);
1819 return iterator1 (*this, rank, i, j, itv, it);
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);
1830 return iterator1 (*this, rank, i, j, itv, it);
1836 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1837 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 ());
1840 vector_const_subiterator_type itv (data ().lower_bound (layout_type::index_M (i, j)));
1841 vector_const_subiterator_type itv_end (data ().end ());
1843 return const_iterator2 (*this, rank, i, j, itv_end, (*(-- itv)).second.end ());
1845 const_subiterator_type it ((*itv).second.lower_bound (layout_type::index_m (i, j)));
1846 const_subiterator_type it_end ((*itv).second.end ());
1848 // advance to the first available major index
1849 size_type M = itv->first;
1854 m = layout_type::size_m(size1_, size2_);
1856 size_type first_j = layout_type::index_m(M,m);
1857 return const_iterator2 (*this, rank, i, first_j, itv, it);
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 ()) {
1864 return const_iterator2 (*this, rank, i, j, itv, it);
1868 return const_iterator2 (*this, rank, i, j, itv, it);
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);
1879 return const_iterator2 (*this, rank, i, j, itv, it);
1885 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
1886 iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) {
1887 BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ());
1889 vector_subiterator_type itv (data ().lower_bound (layout_type::index_M (i, j)));
1890 vector_subiterator_type itv_end (data ().end ());
1892 return iterator2 (*this, rank, i, j, itv_end, (*(-- itv)).second.end ());
1894 subiterator_type it ((*itv).second.lower_bound (layout_type::index_m (i, j)));
1895 subiterator_type it_end ((*itv).second.end ());
1897 // advance to the first available major index
1898 size_type M = itv->first;
1903 m = layout_type::size_m(size1_, size2_);
1905 size_type first_j = layout_type::index_m(M,m);
1906 return iterator2 (*this, rank, i, first_j, itv, it);
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 ()) {
1913 return iterator2 (*this, rank, i, j, itv, it);
1917 return iterator2 (*this, rank, i, j, itv, it);
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);
1928 return iterator2 (*this, rank, i, j, itv, it);
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> {
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;
1945 typedef const_iterator2 dual_iterator_type;
1946 typedef const_reverse_iterator2 dual_reverse_iterator_type;
1948 // Construction and destruction
1951 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), 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) {}
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_) {}
1961 const_iterator1 &operator ++ () {
1962 if (rank_ == 1 && layout_type::fast_i ())
1965 const self_type &m = (*this) ();
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);
1983 const_iterator1 &operator -- () {
1984 if (rank_ == 1 && layout_type::fast_i ())
1987 const self_type &m = (*this) ();
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);
2008 const_reference operator * () const {
2009 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
2010 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
2012 return (*it_).second;
2014 return (*this) () (i_, j_);
2018 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
2020 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2021 typename self_type::
2023 const_iterator2 begin () const {
2024 const self_type &m = (*this) ();
2025 return m.find2 (1, index1 (), 0);
2028 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2029 typename self_type::
2031 const_iterator2 cbegin () const {
2035 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2036 typename self_type::
2038 const_iterator2 end () const {
2039 const self_type &m = (*this) ();
2040 return m.find2 (1, index1 (), m.size2 ());
2043 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2044 typename self_type::
2046 const_iterator2 cend () const {
2050 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2051 typename self_type::
2053 const_reverse_iterator2 rbegin () const {
2054 return const_reverse_iterator2 (end ());
2057 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2058 typename self_type::
2060 const_reverse_iterator2 crbegin () const {
2064 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2065 typename self_type::
2067 const_reverse_iterator2 rend () const {
2068 return const_reverse_iterator2 (begin ());
2071 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2072 typename self_type::
2074 const_reverse_iterator2 crend () const {
2081 size_type index1 () const {
2082 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
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);
2091 size_type index2 () const {
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);
2102 const_iterator1 &operator = (const const_iterator1 &it) {
2103 container_const_reference<self_type>::assign (&it ());
2114 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_;
2120 return i_ == it.i_ && j_ == it.j_;
2128 vector_const_subiterator_type itv_;
2129 const_subiterator_type it_;
2133 const_iterator1 begin1 () const {
2134 return find1 (0, 0, 0);
2137 const_iterator1 cbegin1 () const {
2141 const_iterator1 end1 () const {
2142 return find1 (0, size1_, 0);
2145 const_iterator1 cend1 () const {
2150 public container_reference<mapped_vector_of_mapped_vector>,
2151 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
2152 iterator1, value_type> {
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;
2159 typedef iterator2 dual_iterator_type;
2160 typedef reverse_iterator2 dual_reverse_iterator_type;
2162 // Construction and destruction
2165 container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), 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) {}
2172 iterator1 &operator ++ () {
2173 if (rank_ == 1 && layout_type::fast_i ())
2176 self_type &m = (*this) ();
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);
2194 iterator1 &operator -- () {
2195 if (rank_ == 1 && layout_type::fast_i ())
2198 self_type &m = (*this) ();
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);
2219 reference operator * () const {
2220 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
2221 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
2223 return (*it_).second;
2225 return (*this) ().at_element (i_, j_);
2229 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
2231 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2232 typename self_type::
2234 iterator2 begin () const {
2235 self_type &m = (*this) ();
2236 return m.find2 (1, index1 (), 0);
2239 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2240 typename self_type::
2242 iterator2 end () const {
2243 self_type &m = (*this) ();
2244 return m.find2 (1, index1 (), m.size2 ());
2247 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2248 typename self_type::
2250 reverse_iterator2 rbegin () const {
2251 return reverse_iterator2 (end ());
2254 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2255 typename self_type::
2257 reverse_iterator2 rend () const {
2258 return reverse_iterator2 (begin ());
2264 size_type index1 () const {
2265 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
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);
2274 size_type index2 () const {
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);
2285 iterator1 &operator = (const iterator1 &it) {
2286 container_reference<self_type>::assign (&it ());
2297 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_;
2303 return i_ == it.i_ && j_ == it.j_;
2311 vector_subiterator_type itv_;
2312 subiterator_type it_;
2314 friend class const_iterator1;
2318 iterator1 begin1 () {
2319 return find1 (0, 0, 0);
2323 return find1 (0, size1_, 0);
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> {
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;
2336 typedef const_iterator1 dual_iterator_type;
2337 typedef const_reverse_iterator1 dual_reverse_iterator_type;
2339 // Construction and destruction
2342 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), 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) {}
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_) {}
2352 const_iterator2 &operator ++ () {
2353 if (rank_ == 1 && layout_type::fast_j ())
2356 const self_type &m = (*this) ();
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);
2374 const_iterator2 &operator -- () {
2375 if (rank_ == 1 && layout_type::fast_j ())
2378 const self_type &m = (*this) ();
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);
2399 const_reference operator * () const {
2400 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
2401 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
2403 return (*it_).second;
2405 return (*this) () (i_, j_);
2409 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
2411 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2412 typename self_type::
2414 const_iterator1 begin () const {
2415 const self_type &m = (*this) ();
2416 return m.find1 (1, 0, index2 ());
2419 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2420 typename self_type::
2422 const_iterator1 cbegin () const {
2426 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2427 typename self_type::
2429 const_iterator1 end () const {
2430 const self_type &m = (*this) ();
2431 return m.find1 (1, m.size1 (), index2 ());
2434 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2435 typename self_type::
2437 const_iterator1 cend () const {
2441 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2442 typename self_type::
2444 const_reverse_iterator1 rbegin () const {
2445 return const_reverse_iterator1 (end ());
2448 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2449 typename self_type::
2451 const_reverse_iterator1 crbegin () const {
2455 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2456 typename self_type::
2458 const_reverse_iterator1 rend () const {
2459 return const_reverse_iterator1 (begin ());
2462 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2463 typename self_type::
2465 const_reverse_iterator1 crend () const {
2472 size_type index1 () const {
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);
2481 size_type index2 () const {
2482 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
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);
2493 const_iterator2 &operator = (const const_iterator2 &it) {
2494 container_const_reference<self_type>::assign (&it ());
2505 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_;
2511 return i_ == it.i_ && j_ == it.j_;
2519 vector_const_subiterator_type itv_;
2520 const_subiterator_type it_;
2524 const_iterator2 begin2 () const {
2525 return find2 (0, 0, 0);
2528 const_iterator2 cbegin2 () const {
2532 const_iterator2 end2 () const {
2533 return find2 (0, 0, size2_);
2536 const_iterator2 cend2 () const {
2541 public container_reference<mapped_vector_of_mapped_vector>,
2542 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
2543 iterator2, value_type> {
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;
2550 typedef iterator1 dual_iterator_type;
2551 typedef reverse_iterator1 dual_reverse_iterator_type;
2553 // Construction and destruction
2556 container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), 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) {}
2563 iterator2 &operator ++ () {
2564 if (rank_ == 1 && layout_type::fast_j ())
2567 self_type &m = (*this) ();
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);
2585 iterator2 &operator -- () {
2586 if (rank_ == 1 && layout_type::fast_j ())
2589 self_type &m = (*this) ();
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);
2610 reference operator * () const {
2611 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
2612 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
2614 return (*it_).second;
2616 return (*this) ().at_element (i_, j_);
2620 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
2622 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2623 typename self_type::
2625 iterator1 begin () const {
2626 self_type &m = (*this) ();
2627 return m.find1 (1, 0, index2 ());
2630 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2631 typename self_type::
2633 iterator1 end () const {
2634 self_type &m = (*this) ();
2635 return m.find1 (1, m.size1 (), index2 ());
2638 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2639 typename self_type::
2641 reverse_iterator1 rbegin () const {
2642 return reverse_iterator1 (end ());
2645 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
2646 typename self_type::
2648 reverse_iterator1 rend () const {
2649 return reverse_iterator1 (begin ());
2655 size_type index1 () const {
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);
2664 size_type index2 () const {
2665 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
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);
2676 iterator2 &operator = (const iterator2 &it) {
2677 container_reference<self_type>::assign (&it ());
2688 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_;
2694 return i_ == it.i_ && j_ == it.j_;
2702 vector_subiterator_type itv_;
2703 subiterator_type it_;
2705 friend class const_iterator2;
2709 iterator2 begin2 () {
2710 return find2 (0, 0, 0);
2714 return find2 (0, 0, size2_);
2717 // Reverse iterators
2720 const_reverse_iterator1 rbegin1 () const {
2721 return const_reverse_iterator1 (end1 ());
2724 const_reverse_iterator1 crbegin1 () const {
2728 const_reverse_iterator1 rend1 () const {
2729 return const_reverse_iterator1 (begin1 ());
2732 const_reverse_iterator1 crend1 () const {
2737 reverse_iterator1 rbegin1 () {
2738 return reverse_iterator1 (end1 ());
2741 reverse_iterator1 rend1 () {
2742 return reverse_iterator1 (begin1 ());
2746 const_reverse_iterator2 rbegin2 () const {
2747 return const_reverse_iterator2 (end2 ());
2750 const_reverse_iterator2 crbegin2 () const {
2754 const_reverse_iterator2 rend2 () const {
2755 return const_reverse_iterator2 (begin2 ());
2758 const_reverse_iterator2 crend2 () const {
2763 reverse_iterator2 rbegin2 () {
2764 return reverse_iterator2 (end2 ());
2767 reverse_iterator2 rend2 () {
2768 return reverse_iterator2 (begin2 ());
2772 template<class Archive>
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) {
2782 ar & serialization::make_nvp("data", data_);
2789 static const value_type zero_;
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*/();
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> > {
2802 typedef T &true_reference;
2804 typedef const T *const_pointer;
2805 typedef L layout_type;
2806 typedef compressed_matrix<T, L, IB, IA, TA> self_type;
2808 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
2809 using matrix_container<self_type>::operator ();
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;
2823 typedef sparse_matrix_element<self_type> reference;
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;
2834 // Construction and destruction
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 ();
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 ();
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 ();
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)
2869 reserve(m.nnz(), false);
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;
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 );
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 ();
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);
2901 size_type size1 () const {
2905 size_type size2 () const {
2909 size_type nnz_capacity () const {
2913 size_type nnz () const {
2917 // Storage accessors
2919 static size_type index_base () {
2923 array_size_type filled1 () const {
2927 array_size_type filled2 () const {
2931 const index_array_type &index1_data () const {
2932 return index1_data_;
2935 const index_array_type &index2_data () const {
2936 return index2_data_;
2939 const value_array_type &value_data () const {
2943 void set_filled (const array_size_type& filled1, const array_size_type& filled2) {
2946 storage_invariants ();
2949 index_array_type &index1_data () {
2950 return index1_data_;
2953 index_array_type &index2_data () {
2954 return index2_data_;
2957 value_array_type &value_data () {
2961 void complete_index1_data () {
2962 while (filled1_ <= layout_type::size_M (size1_, size2_)) {
2963 this->index1_data_ [filled1_] = k_based (filled2_);
2971 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_;
2981 void resize (size_type size1, size_type size2, bool preserve = true) {
2982 // FIXME preserve unimplemented
2983 BOOST_UBLAS_CHECK (!preserve, internal_logic ());
2986 capacity_ = restrict_capacity (capacity_);
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 ();
2998 void reserve (size_type non_zeros, bool preserve = true) {
2999 capacity_ = restrict_capacity (non_zeros);
3001 index2_data_.resize (capacity_, size_type ());
3002 value_data_.resize (capacity_, value_type ());
3003 filled2_ = (std::min) (capacity_, filled2_);
3006 index2_data_.resize (capacity_);
3007 value_data_.resize (capacity_);
3010 index1_data_ [filled1_ - 1] = k_based (filled2_);
3012 storage_invariants ();
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));
3021 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)
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))
3032 return &value_data_ [it - index2_data_.begin ()];
3037 const_reference operator () (size_type i, size_type j) const {
3038 const_pointer p = find_element (i, 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);
3055 return insert_element (i, j, value_type/*zero*/());
3057 return reference (*this, i, j);
3061 // Element assignment
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_);
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
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_);
3087 while (element1 + 1 < filled1_) {
3088 ++ index1_data_ [element1 + 1];
3091 storage_invariants ();
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_)
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);
3110 while (index1_data_ [filled1_ - 2] > k_based (filled2_)) {
3111 index1_data_ [filled1_ - 1] = 0;
3114 while (element1 + 1 < filled1_) {
3115 -- index1_data_ [element1 + 1];
3119 storage_invariants ();
3127 index1_data_ [filled1_ - 1] = k_based (filled2_);
3128 storage_invariants ();
3133 compressed_matrix &operator = (const compressed_matrix &m) {
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_;
3144 storage_invariants ();
3147 template<class C> // Container assignment without temporary
3149 compressed_matrix &operator = (const matrix_container<C> &m) {
3150 resize (m ().size1 (), m ().size2 (), false);
3155 compressed_matrix &assign_temporary (compressed_matrix &m) {
3161 compressed_matrix &operator = (const matrix_expression<AE> &ae) {
3162 self_type temporary (ae, capacity_);
3163 return assign_temporary (temporary);
3167 compressed_matrix &assign (const matrix_expression<AE> &ae) {
3168 matrix_assign<scalar_assign> (*this, ae);
3173 compressed_matrix& operator += (const matrix_expression<AE> &ae) {
3174 self_type temporary (*this + ae, capacity_);
3175 return assign_temporary (temporary);
3177 template<class C> // Container assignment without temporary
3179 compressed_matrix &operator += (const matrix_container<C> &m) {
3185 compressed_matrix &plus_assign (const matrix_expression<AE> &ae) {
3186 matrix_assign<scalar_plus_assign> (*this, ae);
3191 compressed_matrix& operator -= (const matrix_expression<AE> &ae) {
3192 self_type temporary (*this - ae, capacity_);
3193 return assign_temporary (temporary);
3195 template<class C> // Container assignment without temporary
3197 compressed_matrix &operator -= (const matrix_container<C> &m) {
3203 compressed_matrix &minus_assign (const matrix_expression<AE> &ae) {
3204 matrix_assign<scalar_minus_assign> (*this, ae);
3209 compressed_matrix& operator *= (const AT &at) {
3210 matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
3215 compressed_matrix& operator /= (const AT &at) {
3216 matrix_assign_scalar<scalar_divides_assign> (*this, at);
3222 void swap (compressed_matrix &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_);
3233 storage_invariants ();
3236 friend void swap (compressed_matrix &m1, compressed_matrix &m2) {
3240 // Back element insertion and erasure
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_);
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 ());
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 ();
3264 BOOST_UBLAS_CHECK (filled1_ > 0 && filled2_ > 0, external_logic ());
3266 while (index1_data_ [filled1_ - 2] > k_based (filled2_)) {
3267 index1_data_ [filled1_ - 1] = 0;
3270 -- index1_data_ [filled1_ - 1];
3271 storage_invariants ();
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;
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 ());
3290 class const_iterator1;
3292 class const_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;
3300 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
3301 const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const {
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_);
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)));
3312 const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ()));
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 ()) {
3320 return const_iterator1 (*this, rank, i, j, itv, it);
3321 i = zero_based (*it);
3324 return const_iterator1 (*this, rank, i, j, itv, it);
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));
3334 return const_iterator1 (*this, rank, i, j, itv, it);
3340 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
3341 iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) {
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_);
3349 subiterator_type it_begin (index2_data_.begin () + zero_based (*itv));
3350 subiterator_type it_end (index2_data_.begin () + zero_based (*(itv + 1)));
3352 subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ()));
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 ()) {
3360 return iterator1 (*this, rank, i, j, itv, it);
3361 i = zero_based (*it);
3364 return iterator1 (*this, rank, i, j, itv, it);
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));
3374 return iterator1 (*this, rank, i, j, itv, it);
3380 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
3381 const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const {
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_);
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)));
3392 const_subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ()));
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 ()) {
3400 return const_iterator2 (*this, rank, i, j, itv, it);
3401 j = zero_based (*it);
3404 return const_iterator2 (*this, rank, i, j, itv, it);
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));
3414 return const_iterator2 (*this, rank, i, j, itv, it);
3420 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
3421 iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) {
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_);
3429 subiterator_type it_begin (index2_data_.begin () + zero_based (*itv));
3430 subiterator_type it_end (index2_data_.begin () + zero_based (*(itv + 1)));
3432 subiterator_type it (detail::lower_bound (it_begin, it_end, k_based (address2), std::less<size_type> ()));
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 ()) {
3440 return iterator2 (*this, rank, i, j, itv, it);
3441 j = zero_based (*it);
3444 return iterator2 (*this, rank, i, j, itv, it);
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));
3454 return iterator2 (*this, rank, i, j, itv, it);
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> {
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;
3472 typedef const_iterator2 dual_iterator_type;
3473 typedef const_reverse_iterator2 dual_reverse_iterator_type;
3475 // Construction and destruction
3478 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), 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) {}
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_) {}
3488 const_iterator1 &operator ++ () {
3489 if (rank_ == 1 && layout_type::fast_i ())
3494 *this = (*this) ().find1 (rank_, i_, j_, 1);
3499 const_iterator1 &operator -- () {
3500 if (rank_ == 1 && layout_type::fast_i ())
3505 *this = (*this) ().find1 (rank_, i_, j_, -1);
3512 const_reference operator * () const {
3513 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
3514 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
3516 return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
3518 return (*this) () (i_, j_);
3522 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
3524 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3525 typename self_type::
3527 const_iterator2 begin () const {
3528 const self_type &m = (*this) ();
3529 return m.find2 (1, index1 (), 0);
3532 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3533 typename self_type::
3535 const_iterator2 cbegin () const {
3539 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3540 typename self_type::
3542 const_iterator2 end () const {
3543 const self_type &m = (*this) ();
3544 return m.find2 (1, index1 (), m.size2 ());
3547 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3548 typename self_type::
3550 const_iterator2 cend () const {
3554 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3555 typename self_type::
3557 const_reverse_iterator2 rbegin () const {
3558 return const_reverse_iterator2 (end ());
3561 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3562 typename self_type::
3564 const_reverse_iterator2 crbegin () const {
3568 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3569 typename self_type::
3571 const_reverse_iterator2 rend () const {
3572 return const_reverse_iterator2 (begin ());
3575 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3576 typename self_type::
3578 const_reverse_iterator2 crend () const {
3585 size_type index1 () const {
3586 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
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_));
3595 size_type index2 () const {
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_));
3606 const_iterator1 &operator = (const const_iterator1 &it) {
3607 container_const_reference<self_type>::assign (&it ());
3618 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_;
3624 return i_ == it.i_ && j_ == it.j_;
3632 vector_const_subiterator_type itv_;
3633 const_subiterator_type it_;
3637 const_iterator1 begin1 () const {
3638 return find1 (0, 0, 0);
3641 const_iterator1 cbegin1 () const {
3645 const_iterator1 end1 () const {
3646 return find1 (0, size1_, 0);
3649 const_iterator1 cend1 () const {
3654 public container_reference<compressed_matrix>,
3655 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
3656 iterator1, value_type> {
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;
3663 typedef iterator2 dual_iterator_type;
3664 typedef reverse_iterator2 dual_reverse_iterator_type;
3666 // Construction and destruction
3669 container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), 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) {}
3676 iterator1 &operator ++ () {
3677 if (rank_ == 1 && layout_type::fast_i ())
3682 *this = (*this) ().find1 (rank_, i_, j_, 1);
3687 iterator1 &operator -- () {
3688 if (rank_ == 1 && layout_type::fast_i ())
3693 *this = (*this) ().find1 (rank_, i_, j_, -1);
3700 reference operator * () const {
3701 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
3702 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
3704 return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
3706 return (*this) ().at_element (i_, j_);
3710 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
3712 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3713 typename self_type::
3715 iterator2 begin () const {
3716 self_type &m = (*this) ();
3717 return m.find2 (1, index1 (), 0);
3720 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3721 typename self_type::
3723 iterator2 end () const {
3724 self_type &m = (*this) ();
3725 return m.find2 (1, index1 (), m.size2 ());
3728 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3729 typename self_type::
3731 reverse_iterator2 rbegin () const {
3732 return reverse_iterator2 (end ());
3735 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3736 typename self_type::
3738 reverse_iterator2 rend () const {
3739 return reverse_iterator2 (begin ());
3745 size_type index1 () const {
3746 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
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_));
3755 size_type index2 () const {
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_));
3766 iterator1 &operator = (const iterator1 &it) {
3767 container_reference<self_type>::assign (&it ());
3778 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_;
3784 return i_ == it.i_ && j_ == it.j_;
3792 vector_subiterator_type itv_;
3793 subiterator_type it_;
3795 friend class const_iterator1;
3799 iterator1 begin1 () {
3800 return find1 (0, 0, 0);
3804 return find1 (0, size1_, 0);
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> {
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;
3817 typedef const_iterator1 dual_iterator_type;
3818 typedef const_reverse_iterator1 dual_reverse_iterator_type;
3820 // Construction and destruction
3823 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), 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) {}
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_) {}
3833 const_iterator2 &operator ++ () {
3834 if (rank_ == 1 && layout_type::fast_j ())
3839 *this = (*this) ().find2 (rank_, i_, j_, 1);
3844 const_iterator2 &operator -- () {
3845 if (rank_ == 1 && layout_type::fast_j ())
3850 *this = (*this) ().find2 (rank_, i_, j_, -1);
3857 const_reference operator * () const {
3858 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
3859 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
3861 return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
3863 return (*this) () (i_, j_);
3867 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
3869 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3870 typename self_type::
3872 const_iterator1 begin () const {
3873 const self_type &m = (*this) ();
3874 return m.find1 (1, 0, index2 ());
3877 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3878 typename self_type::
3880 const_iterator1 cbegin () const {
3884 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3885 typename self_type::
3887 const_iterator1 end () const {
3888 const self_type &m = (*this) ();
3889 return m.find1 (1, m.size1 (), index2 ());
3892 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3893 typename self_type::
3895 const_iterator1 cend () const {
3899 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3900 typename self_type::
3902 const_reverse_iterator1 rbegin () const {
3903 return const_reverse_iterator1 (end ());
3906 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3907 typename self_type::
3909 const_reverse_iterator1 crbegin () const {
3913 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3914 typename self_type::
3916 const_reverse_iterator1 rend () const {
3917 return const_reverse_iterator1 (begin ());
3920 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
3921 typename self_type::
3923 const_reverse_iterator1 crend () const {
3930 size_type index1 () const {
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_));
3939 size_type index2 () const {
3940 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
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_));
3951 const_iterator2 &operator = (const const_iterator2 &it) {
3952 container_const_reference<self_type>::assign (&it ());
3963 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_;
3969 return i_ == it.i_ && j_ == it.j_;
3977 vector_const_subiterator_type itv_;
3978 const_subiterator_type it_;
3982 const_iterator2 begin2 () const {
3983 return find2 (0, 0, 0);
3986 const_iterator2 cbegin2 () const {
3990 const_iterator2 end2 () const {
3991 return find2 (0, 0, size2_);
3994 const_iterator2 cend2 () const {
3999 public container_reference<compressed_matrix>,
4000 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
4001 iterator2, value_type> {
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;
4008 typedef iterator1 dual_iterator_type;
4009 typedef reverse_iterator1 dual_reverse_iterator_type;
4011 // Construction and destruction
4014 container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), 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) {}
4021 iterator2 &operator ++ () {
4022 if (rank_ == 1 && layout_type::fast_j ())
4027 *this = (*this) ().find2 (rank_, i_, j_, 1);
4032 iterator2 &operator -- () {
4033 if (rank_ == 1 && layout_type::fast_j ())
4038 *this = (*this) ().find2 (rank_, i_, j_, -1);
4045 reference operator * () const {
4046 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
4047 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
4049 return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
4051 return (*this) ().at_element (i_, j_);
4055 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
4057 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
4058 typename self_type::
4060 iterator1 begin () const {
4061 self_type &m = (*this) ();
4062 return m.find1 (1, 0, index2 ());
4065 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
4066 typename self_type::
4068 iterator1 end () const {
4069 self_type &m = (*this) ();
4070 return m.find1 (1, m.size1 (), index2 ());
4073 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
4074 typename self_type::
4076 reverse_iterator1 rbegin () const {
4077 return reverse_iterator1 (end ());
4080 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
4081 typename self_type::
4083 reverse_iterator1 rend () const {
4084 return reverse_iterator1 (begin ());
4090 size_type index1 () const {
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_));
4099 size_type index2 () const {
4100 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
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_));
4111 iterator2 &operator = (const iterator2 &it) {
4112 container_reference<self_type>::assign (&it ());
4123 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_;
4129 return i_ == it.i_ && j_ == it.j_;
4137 vector_subiterator_type itv_;
4138 subiterator_type it_;
4140 friend class const_iterator2;
4144 iterator2 begin2 () {
4145 return find2 (0, 0, 0);
4149 return find2 (0, 0, size2_);
4152 // Reverse iterators
4155 const_reverse_iterator1 rbegin1 () const {
4156 return const_reverse_iterator1 (end1 ());
4159 const_reverse_iterator1 crbegin1 () const {
4163 const_reverse_iterator1 rend1 () const {
4164 return const_reverse_iterator1 (begin1 ());
4167 const_reverse_iterator1 crend1 () const {
4172 reverse_iterator1 rbegin1 () {
4173 return reverse_iterator1 (end1 ());
4176 reverse_iterator1 rend1 () {
4177 return reverse_iterator1 (begin1 ());
4181 const_reverse_iterator2 rbegin2 () const {
4182 return const_reverse_iterator2 (end2 ());
4185 const_reverse_iterator2 crbegin2 () const {
4189 const_reverse_iterator2 rend2 () const {
4190 return const_reverse_iterator2 (begin2 ());
4193 const_reverse_iterator2 crend2 () const {
4198 reverse_iterator2 rbegin2 () {
4199 return reverse_iterator2 (end2 ());
4202 reverse_iterator2 rend2 () {
4203 return reverse_iterator2 (begin2 ());
4207 template<class Archive>
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) {
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();
4227 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 ());
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_;
4247 static size_type zero_based (size_type k_based_index) {
4248 return k_based_index - IB;
4251 static size_type k_based (size_type zero_based_index) {
4252 return zero_based_index + IB;
4255 friend class iterator1;
4256 friend class iterator2;
4257 friend class const_iterator1;
4258 friend class const_iterator2;
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*/();
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> > {
4271 typedef T &true_reference;
4273 typedef const T *const_pointer;
4274 typedef L layout_type;
4275 typedef coordinate_matrix<T, L, IB, IA, TA> self_type;
4277 #ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
4278 using matrix_container<self_type>::operator ();
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;
4291 typedef sparse_matrix_element<self_type> reference;
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;
4302 // Construction and destruction
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 ();
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 ();
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 ();
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);
4340 size_type size1 () const {
4344 size_type size2 () const {
4348 size_type nnz_capacity () const {
4352 size_type nnz () const {
4356 // Storage accessors
4358 static size_type index_base () {
4362 array_size_type filled () const {
4366 const index_array_type &index1_data () const {
4367 return index1_data_;
4370 const index_array_type &index2_data () const {
4371 return index2_data_;
4374 const value_array_type &value_data () const {
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;
4383 sorted_ = (sorted_filled_ == filled);
4385 storage_invariants ();
4388 index_array_type &index1_data () {
4389 return index1_data_;
4392 index_array_type &index2_data () {
4393 return index2_data_;
4396 value_array_type &value_data () {
4403 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
4411 void resize (size_type size1, size_type size2, bool preserve = true) {
4412 // FIXME preserve unimplemented
4413 BOOST_UBLAS_CHECK (!preserve, internal_logic ());
4416 capacity_ = restrict_capacity (capacity_);
4417 index1_data_.resize (capacity_);
4418 index2_data_.resize (capacity_);
4419 value_data_.resize (capacity_);
4421 sorted_filled_ = filled_;
4423 storage_invariants ();
4428 void reserve (array_size_type non_zeros, bool preserve = true) {
4429 sort (); // remove duplicate elements
4430 capacity_ = restrict_capacity (non_zeros);
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_);
4438 index1_data_.resize (capacity_);
4439 index2_data_.resize (capacity_);
4440 value_data_.resize (capacity_);
4443 sorted_filled_ = filled_;
4444 storage_invariants ();
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));
4453 const_pointer find_element (size_type i, size_type j) const {
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)
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))
4466 return &value_data_ [it - index2_data_.begin ()];
4471 const_reference operator () (size_type i, size_type j) const {
4472 const_pointer p = find_element (i, j);
4479 reference operator () (size_type i, size_type j) {
4480 #ifndef BOOST_UBLAS_STRICT_MATRIX_SPARSE
4481 pointer p = find_element (i, j);
4485 return insert_element (i, j, value_type/*zero*/());
4487 return reference (*this, i, j);
4491 // Element assignment
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;
4504 storage_invariants ();
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];
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);
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);
4530 sorted_filled_ = filled_;
4532 storage_invariants ();
4539 sorted_filled_ = filled_;
4541 storage_invariants ();
4546 coordinate_matrix &operator = (const coordinate_matrix &m) {
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 ());
4561 storage_invariants ();
4564 template<class C> // Container assignment without temporary
4566 coordinate_matrix &operator = (const matrix_container<C> &m) {
4567 resize (m ().size1 (), m ().size2 (), false);
4572 coordinate_matrix &assign_temporary (coordinate_matrix &m) {
4578 coordinate_matrix &operator = (const matrix_expression<AE> &ae) {
4579 self_type temporary (ae, capacity_);
4580 return assign_temporary (temporary);
4584 coordinate_matrix &assign (const matrix_expression<AE> &ae) {
4585 matrix_assign<scalar_assign> (*this, ae);
4590 coordinate_matrix& operator += (const matrix_expression<AE> &ae) {
4591 self_type temporary (*this + ae, capacity_);
4592 return assign_temporary (temporary);
4594 template<class C> // Container assignment without temporary
4596 coordinate_matrix &operator += (const matrix_container<C> &m) {
4602 coordinate_matrix &plus_assign (const matrix_expression<AE> &ae) {
4603 matrix_assign<scalar_plus_assign> (*this, ae);
4608 coordinate_matrix& operator -= (const matrix_expression<AE> &ae) {
4609 self_type temporary (*this - ae, capacity_);
4610 return assign_temporary (temporary);
4612 template<class C> // Container assignment without temporary
4614 coordinate_matrix &operator -= (const matrix_container<C> &m) {
4620 coordinate_matrix &minus_assign (const matrix_expression<AE> &ae) {
4621 matrix_assign<scalar_minus_assign> (*this, ae);
4626 coordinate_matrix& operator *= (const AT &at) {
4627 matrix_assign_scalar<scalar_multiplies_assign> (*this, at);
4632 coordinate_matrix& operator /= (const AT &at) {
4633 matrix_assign_scalar<scalar_divides_assign> (*this, at);
4639 void swap (coordinate_matrix &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_);
4651 storage_invariants ();
4654 friend void swap (coordinate_matrix &m1, coordinate_matrix &m2) {
4658 // replacement if STL lower bound algorithm for use of inplace_merge
4659 array_size_type lower_bound (array_size_type beg, array_size_type end, array_size_type target) const {
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])))) {
4673 // specialized replacement of STL inplace_merge to avoid compilation
4674 // problems with respect to the array_triple iterator
4675 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;
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])))
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]);
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);
4693 rig_mid = (mid + end) / 2;
4694 lef_mid = lower_bound(beg, mid, rig_mid);
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);
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);
4706 // Sorting and summation of duplicates
4708 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_);
4718 const typename array_triple::iterator iunsorted = ita.begin ();
4719 std::sort (iunsorted, ita.end ());
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]) {
4728 index1_data_ [filled] = index1_data_ [i];
4729 index2_data_ [filled] = index2_data_ [i];
4730 value_data_ [filled] = value_data_ [i];
4733 value_data_ [filled] += value_data_ [i];
4736 filled_ = filled + 1;
4737 sorted_filled_ = filled_;
4739 storage_invariants ();
4743 // Back element insertion and erasure
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_ &&
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;
4761 sorted_filled_ = filled_;
4762 storage_invariants ();
4766 // ISSUE invariants could be simpilfied if sorted required as precondition
4767 BOOST_UBLAS_CHECK (filled_ > 0, external_logic ());
4769 sorted_filled_ = (std::min) (sorted_filled_, filled_);
4770 sorted_ = sorted_filled_ = filled_;
4771 storage_invariants ();
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;
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 ());
4790 class const_iterator1;
4792 class const_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;
4800 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
4801 const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const {
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> ()));
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 ()));
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 ()));
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 ()) {
4821 return const_iterator1 (*this, rank, i, j, itv, it);
4822 i = zero_based (*it);
4825 return const_iterator1 (*this, rank, i, j, itv, it);
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));
4835 return const_iterator1 (*this, rank, i, j, itv, it);
4841 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
4842 iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) {
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> ()));
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 ()));
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 ()));
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 ()) {
4862 return iterator1 (*this, rank, i, j, itv, it);
4863 i = zero_based (*it);
4866 return iterator1 (*this, rank, i, j, itv, it);
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));
4876 return iterator1 (*this, rank, i, j, itv, it);
4882 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
4883 const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const {
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> ()));
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 ()));
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 ()));
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 ()) {
4903 return const_iterator2 (*this, rank, i, j, itv, it);
4904 j = zero_based (*it);
4907 return const_iterator2 (*this, rank, i, j, itv, it);
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));
4917 return const_iterator2 (*this, rank, i, j, itv, it);
4923 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
4924 iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) {
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> ()));
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 ()));
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 ()));
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 ()) {
4944 return iterator2 (*this, rank, i, j, itv, it);
4945 j = zero_based (*it);
4948 return iterator2 (*this, rank, i, j, itv, it);
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));
4958 return iterator2 (*this, rank, i, j, itv, it);
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> {
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;
4976 typedef const_iterator2 dual_iterator_type;
4977 typedef const_reverse_iterator2 dual_reverse_iterator_type;
4979 // Construction and destruction
4982 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), 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) {}
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_) {}
4992 const_iterator1 &operator ++ () {
4993 if (rank_ == 1 && layout_type::fast_i ())
4998 *this = (*this) ().find1 (rank_, i_, j_, 1);
5003 const_iterator1 &operator -- () {
5004 if (rank_ == 1 && layout_type::fast_i ())
5009 *this = (*this) ().find1 (rank_, i_, j_, -1);
5016 const_reference operator * () const {
5017 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
5018 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
5020 return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
5022 return (*this) () (i_, j_);
5026 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
5028 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5029 typename self_type::
5031 const_iterator2 begin () const {
5032 const self_type &m = (*this) ();
5033 return m.find2 (1, index1 (), 0);
5036 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5037 typename self_type::
5039 const_iterator2 cbegin () const {
5043 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5044 typename self_type::
5046 const_iterator2 end () const {
5047 const self_type &m = (*this) ();
5048 return m.find2 (1, index1 (), m.size2 ());
5051 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5052 typename self_type::
5054 const_iterator2 cend () const {
5058 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5059 typename self_type::
5061 const_reverse_iterator2 rbegin () const {
5062 return const_reverse_iterator2 (end ());
5065 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5066 typename self_type::
5068 const_reverse_iterator2 crbegin () const {
5072 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5073 typename self_type::
5075 const_reverse_iterator2 rend () const {
5076 return const_reverse_iterator2 (begin ());
5079 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5080 typename self_type::
5082 const_reverse_iterator2 crend () const {
5089 size_type index1 () const {
5090 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
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_));
5099 size_type index2 () const {
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_));
5110 const_iterator1 &operator = (const const_iterator1 &it) {
5111 container_const_reference<self_type>::assign (&it ());
5122 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_;
5128 return i_ == it.i_ && j_ == it.j_;
5136 vector_const_subiterator_type itv_;
5137 const_subiterator_type it_;
5141 const_iterator1 begin1 () const {
5142 return find1 (0, 0, 0);
5145 const_iterator1 cbegin1 () const {
5149 const_iterator1 end1 () const {
5150 return find1 (0, size1_, 0);
5153 const_iterator1 cend1 () const {
5158 public container_reference<coordinate_matrix>,
5159 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
5160 iterator1, value_type> {
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;
5167 typedef iterator2 dual_iterator_type;
5168 typedef reverse_iterator2 dual_reverse_iterator_type;
5170 // Construction and destruction
5173 container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), 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) {}
5180 iterator1 &operator ++ () {
5181 if (rank_ == 1 && layout_type::fast_i ())
5186 *this = (*this) ().find1 (rank_, i_, j_, 1);
5191 iterator1 &operator -- () {
5192 if (rank_ == 1 && layout_type::fast_i ())
5197 *this = (*this) ().find1 (rank_, i_, j_, -1);
5204 reference operator * () const {
5205 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
5206 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
5208 return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
5210 return (*this) ().at_element (i_, j_);
5214 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
5216 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5217 typename self_type::
5219 iterator2 begin () const {
5220 self_type &m = (*this) ();
5221 return m.find2 (1, index1 (), 0);
5224 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5225 typename self_type::
5227 iterator2 end () const {
5228 self_type &m = (*this) ();
5229 return m.find2 (1, index1 (), m.size2 ());
5232 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5233 typename self_type::
5235 reverse_iterator2 rbegin () const {
5236 return reverse_iterator2 (end ());
5239 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5240 typename self_type::
5242 reverse_iterator2 rend () const {
5243 return reverse_iterator2 (begin ());
5249 size_type index1 () const {
5250 BOOST_UBLAS_CHECK (*this != (*this) ().find1 (0, (*this) ().size1 (), j_), bad_index ());
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_));
5259 size_type index2 () const {
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_));
5270 iterator1 &operator = (const iterator1 &it) {
5271 container_reference<self_type>::assign (&it ());
5282 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_;
5288 return i_ == it.i_ && j_ == it.j_;
5296 vector_subiterator_type itv_;
5297 subiterator_type it_;
5299 friend class const_iterator1;
5303 iterator1 begin1 () {
5304 return find1 (0, 0, 0);
5308 return find1 (0, size1_, 0);
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> {
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;
5321 typedef const_iterator1 dual_iterator_type;
5322 typedef const_reverse_iterator1 dual_reverse_iterator_type;
5324 // Construction and destruction
5327 container_const_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), 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) {}
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_) {}
5337 const_iterator2 &operator ++ () {
5338 if (rank_ == 1 && layout_type::fast_j ())
5343 *this = (*this) ().find2 (rank_, i_, j_, 1);
5348 const_iterator2 &operator -- () {
5349 if (rank_ == 1 && layout_type::fast_j ())
5354 *this = (*this) ().find2 (rank_, i_, j_, -1);
5361 const_reference operator * () const {
5362 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
5363 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
5365 return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
5367 return (*this) () (i_, j_);
5371 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
5373 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5374 typename self_type::
5376 const_iterator1 begin () const {
5377 const self_type &m = (*this) ();
5378 return m.find1 (1, 0, index2 ());
5381 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5382 typename self_type::
5384 const_iterator1 cbegin () const {
5388 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5389 typename self_type::
5391 const_iterator1 end () const {
5392 const self_type &m = (*this) ();
5393 return m.find1 (1, m.size1 (), index2 ());
5396 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5397 typename self_type::
5399 const_iterator1 cend () const {
5403 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5404 typename self_type::
5406 const_reverse_iterator1 rbegin () const {
5407 return const_reverse_iterator1 (end ());
5410 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5411 typename self_type::
5413 const_reverse_iterator1 crbegin () const {
5417 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5418 typename self_type::
5420 const_reverse_iterator1 rend () const {
5421 return const_reverse_iterator1 (begin ());
5424 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5425 typename self_type::
5427 const_reverse_iterator1 crend () const {
5434 size_type index1 () const {
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_));
5443 size_type index2 () const {
5444 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
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_));
5455 const_iterator2 &operator = (const const_iterator2 &it) {
5456 container_const_reference<self_type>::assign (&it ());
5467 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_;
5473 return i_ == it.i_ && j_ == it.j_;
5481 vector_const_subiterator_type itv_;
5482 const_subiterator_type it_;
5486 const_iterator2 begin2 () const {
5487 return find2 (0, 0, 0);
5490 const_iterator2 cbegin2 () const {
5494 const_iterator2 end2 () const {
5495 return find2 (0, 0, size2_);
5498 const_iterator2 cend2 () const {
5503 public container_reference<coordinate_matrix>,
5504 public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
5505 iterator2, value_type> {
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;
5512 typedef iterator1 dual_iterator_type;
5513 typedef reverse_iterator1 dual_reverse_iterator_type;
5515 // Construction and destruction
5518 container_reference<self_type> (), rank_ (), i_ (), j_ (), itv_ (), 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) {}
5525 iterator2 &operator ++ () {
5526 if (rank_ == 1 && layout_type::fast_j ())
5531 *this = (*this) ().find2 (rank_, i_, j_, 1);
5536 iterator2 &operator -- () {
5537 if (rank_ == 1 && layout_type::fast_j ())
5542 *this = (*this) ().find2 (rank_, i_, j_, -1);
5549 reference operator * () const {
5550 BOOST_UBLAS_CHECK (index1 () < (*this) ().size1 (), bad_index ());
5551 BOOST_UBLAS_CHECK (index2 () < (*this) ().size2 (), bad_index ());
5553 return (*this) ().value_data_ [it_ - (*this) ().index2_data_.begin ()];
5555 return (*this) ().at_element (i_, j_);
5559 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
5561 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5562 typename self_type::
5564 iterator1 begin () const {
5565 self_type &m = (*this) ();
5566 return m.find1 (1, 0, index2 ());
5569 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5570 typename self_type::
5572 iterator1 end () const {
5573 self_type &m = (*this) ();
5574 return m.find1 (1, m.size1 (), index2 ());
5577 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5578 typename self_type::
5580 reverse_iterator1 rbegin () const {
5581 return reverse_iterator1 (end ());
5584 #ifdef BOOST_UBLAS_MSVC_NESTED_CLASS_RELATION
5585 typename self_type::
5587 reverse_iterator1 rend () const {
5588 return reverse_iterator1 (begin ());
5594 size_type index1 () const {
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_));
5603 size_type index2 () const {
5604 BOOST_UBLAS_CHECK (*this != (*this) ().find2 (0, i_, (*this) ().size2 ()), bad_index ());
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_));
5615 iterator2 &operator = (const iterator2 &it) {
5616 container_reference<self_type>::assign (&it ());
5627 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_;
5633 return i_ == it.i_ && j_ == it.j_;
5641 vector_subiterator_type itv_;
5642 subiterator_type it_;
5644 friend class const_iterator2;
5648 iterator2 begin2 () {
5649 return find2 (0, 0, 0);
5653 return find2 (0, 0, size2_);
5656 // Reverse iterators
5659 const_reverse_iterator1 rbegin1 () const {
5660 return const_reverse_iterator1 (end1 ());
5663 const_reverse_iterator1 crbegin1 () const {
5667 const_reverse_iterator1 rend1 () const {
5668 return const_reverse_iterator1 (begin1 ());
5671 const_reverse_iterator1 crend1 () const {
5676 reverse_iterator1 rbegin1 () {
5677 return reverse_iterator1 (end1 ());
5680 reverse_iterator1 rend1 () {
5681 return reverse_iterator1 (begin1 ());
5685 const_reverse_iterator2 rbegin2 () const {
5686 return const_reverse_iterator2 (end2 ());
5689 const_reverse_iterator2 crbegin2 () const {
5693 const_reverse_iterator2 rend2 () const {
5694 return const_reverse_iterator2 (begin2 ());
5697 const_reverse_iterator2 crend2 () const {
5702 reverse_iterator2 rbegin2 () {
5703 return reverse_iterator2 (end2 ());
5706 reverse_iterator2 rend2 () {
5707 return reverse_iterator2 (begin2 ());
5711 template<class Archive>
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) {
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();
5732 void storage_invariants () const
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 ());
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_;
5754 static size_type zero_based (size_type k_based_index) {
5755 return k_based_index - IB;
5758 static size_type k_based (size_type zero_based_index) {
5759 return zero_based_index + IB;
5762 friend class iterator1;
5763 friend class iterator2;
5764 friend class const_iterator1;
5765 friend class const_iterator2;
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*/();