From bb092c5e461ff566d67de0adc0d6b83ae2cf6c67 Mon Sep 17 00:00:00 2001 From: Michael Stevens Date: Wed, 20 Oct 2004 07:46:32 +0000 Subject: [PATCH] FIX nasty coordinate_matrix insert bug Back to sumation of duplicates [SVN r25793] --- include/boost/numeric/ublas/matrix_sparse.hpp | 78 +++++++++++-------- include/boost/numeric/ublas/vector_sparse.hpp | 49 ++++++------ 2 files changed, 73 insertions(+), 54 deletions(-) diff --git a/include/boost/numeric/ublas/matrix_sparse.hpp b/include/boost/numeric/ublas/matrix_sparse.hpp index fd0b9cd3..4ed8eda1 100644 --- a/include/boost/numeric/ublas/matrix_sparse.hpp +++ b/include/boost/numeric/ublas/matrix_sparse.hpp @@ -270,7 +270,7 @@ namespace boost { namespace numeric { namespace ublas { sparse_matrix (size_type size1, size_type size2, size_type non_zeros = 0): matrix_expression (), size1_ (size1), size2_ (size2), data_ () { - detail::map_reserve (data (), max_nz (non_zeros)); + detail::map_reserve (data (), restrict_nz (non_zeros)); } BOOST_UBLAS_INLINE sparse_matrix (const sparse_matrix &m): @@ -281,7 +281,7 @@ namespace boost { namespace numeric { namespace ublas { sparse_matrix (const matrix_expression &ae, size_type non_zeros = 0): matrix_expression (), size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), data_ () { - detail::map_reserve (data (), max_nz (non_zeros)); + detail::map_reserve (data (), restrict_nz (non_zeros)); matrix_assign (scalar_assign (), *this, ae); } @@ -310,7 +310,7 @@ namespace boost { namespace numeric { namespace ublas { // Resizing private: BOOST_UBLAS_INLINE - size_type max_nz (size_type non_zeros) const { + size_type restrict_nz (size_type non_zeros) const { non_zeros = (std::max) (non_zeros, (std::min) (size1_, size2_)); // Guarding against overflow. // Thanks to Alexei Novakov for the hint. @@ -332,7 +332,7 @@ namespace boost { namespace numeric { namespace ublas { // Reserving BOOST_UBLAS_INLINE void reserve (size_type non_zeros, bool preserve = true) { - detail::map_reserve (data (), max_nz (non_zeros)); + detail::map_reserve (data (), restrict_nz (non_zeros)); } // Proxy support @@ -2495,7 +2495,7 @@ namespace boost { namespace numeric { namespace ublas { BOOST_UBLAS_INLINE compressed_matrix (): matrix_expression (), - size1_ (0), size2_ (0), non_zeros_ (max_nz (0)), + size1_ (0), size2_ (0), non_zeros_ (restrict_nz (0)), filled1_ (1), filled2_ (0), index1_data_ (functor_type::size1 (size1_, size2_) + 1), index2_data_ (non_zeros_), value_data_ (non_zeros_) { @@ -2504,7 +2504,7 @@ namespace boost { namespace numeric { namespace ublas { BOOST_UBLAS_INLINE compressed_matrix (size_type size1, size_type size2, size_type non_zeros = 0): matrix_expression (), - size1_ (size1), size2_ (size2), non_zeros_ (max_nz (non_zeros)), + size1_ (size1), size2_ (size2), non_zeros_ (restrict_nz (non_zeros)), filled1_ (1), filled2_ (0), index1_data_ (functor_type::size1 (size1_, size2_) + 1), index2_data_ (non_zeros_), value_data_ (non_zeros_) { @@ -2523,7 +2523,7 @@ namespace boost { namespace numeric { namespace ublas { BOOST_UBLAS_INLINE compressed_matrix (const matrix_expression &ae, size_type non_zeros = 0): matrix_expression (), - size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), non_zeros_ (max_nz (non_zeros)), + size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), non_zeros_ (restrict_nz (non_zeros)), filled1_ (1), filled2_ (0), index1_data_ (functor_type::size1 (ae ().size1 (), ae ().size2 ()) + 1), index2_data_ (non_zeros_), value_data_ (non_zeros_) { @@ -2588,7 +2588,7 @@ namespace boost { namespace numeric { namespace ublas { // Resizing private: BOOST_UBLAS_INLINE - size_type max_nz (size_type non_zeros) const { + size_type restrict_nz (size_type non_zeros) const { non_zeros = (std::max) (non_zeros, (std::min) (size1_, size2_)); // Guarding against overflow. // Thanks to Alexei Novakov for the hint. @@ -2604,7 +2604,7 @@ namespace boost { namespace numeric { namespace ublas { BOOST_UBLAS_CHECK (!preserve, internal_logic ()); size1_ = size1; size2_ = size2; - non_zeros_ = max_nz (non_zeros_); + non_zeros_ = restrict_nz (non_zeros_); filled1_ = 1; filled2_ = 0; index1_data ().resize (functor_type::size1 (size1_, size2_) + 1); @@ -2616,7 +2616,7 @@ namespace boost { namespace numeric { namespace ublas { // Reserving BOOST_UBLAS_INLINE void reserve (size_type non_zeros, bool preserve = true) { - non_zeros_ = max_nz (non_zeros); + non_zeros_ = restrict_nz (non_zeros); if (preserve) { index2_data ().resize (non_zeros_, size_type ()); value_data ().resize (non_zeros_, value_type ()); @@ -3833,14 +3833,14 @@ namespace boost { namespace numeric { namespace ublas { BOOST_UBLAS_INLINE coordinate_matrix (): matrix_expression (), - size1_ (0), size2_ (0), non_zeros_ (max_nz (0)), + size1_ (0), size2_ (0), non_zeros_ (restrict_nz (0)), filled_ (0), sorted_ (true), index1_data_ (non_zeros_), index2_data_ (non_zeros_), value_data_ (non_zeros_) {} BOOST_UBLAS_INLINE coordinate_matrix (size_type size1, size_type size2, size_type non_zeros = 0): matrix_expression (), - size1_ (size1), size2_ (size2), non_zeros_ (max_nz (non_zeros)), + size1_ (size1), size2_ (size2), non_zeros_ (restrict_nz (non_zeros)), filled_ (0), sorted_ (true), index1_data_ (non_zeros_), index2_data_ (non_zeros_), value_data_ (non_zeros_) { @@ -3857,7 +3857,7 @@ namespace boost { namespace numeric { namespace ublas { BOOST_UBLAS_INLINE coordinate_matrix (const matrix_expression &ae, size_type non_zeros = 0): matrix_expression (), - size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), non_zeros_ (max_nz (non_zeros)), + size1_ (ae ().size1 ()), size2_ (ae ().size2 ()), non_zeros_ (restrict_nz (non_zeros)), filled_ (0), sorted_ (true), index1_data_ (non_zeros_), index2_data_ (non_zeros_), value_data_ (non_zeros_) { @@ -3913,13 +3913,10 @@ namespace boost { namespace numeric { namespace ublas { // Resizing private: BOOST_UBLAS_INLINE - size_type max_nz (size_type non_zeros) const { + size_type restrict_nz (size_type non_zeros) const { + // minimum non_zeros non_zeros = (std::max) (non_zeros, (std::min) (size1_, size2_)); - // Guarding against overflow. - // Thanks to Alexei Novakov for the hint. - // non_zeros_ = (std::min) (non_zeros, size1_ * size2_); - if (size1_ > 0 && non_zeros / size1_ >= size2_) - non_zeros = size1_ * size2_; + // ISSUE no maximum as coordinate may contain inserted duplicates return non_zeros; } public: @@ -3929,18 +3926,19 @@ namespace boost { namespace numeric { namespace ublas { BOOST_UBLAS_CHECK (!preserve, internal_logic ()); size1_ = size1; size2_ = size2; - non_zeros_ = max_nz (non_zeros_); + non_zeros_ = restrict_nz (non_zeros_); index1_data ().resize (non_zeros_); index2_data ().resize (non_zeros_); value_data ().resize (non_zeros_); filled_ = 0; + BOOST_UBLAS_CHECK (filled_ <= non_zeros, internal_logic ()); } // Reserving BOOST_UBLAS_INLINE void reserve (size_type non_zeros, bool preserve = true) { sort (); // remove duplicate elements - non_zeros_ = max_nz (non_zeros); + non_zeros_ = restrict_nz (non_zeros); if (preserve) { index1_data ().resize (non_zeros_, size_type ()); index2_data ().resize (non_zeros_, size_type ()); @@ -3953,6 +3951,7 @@ namespace boost { namespace numeric { namespace ublas { value_data ().resize (non_zeros_); filled_ = 0; } + BOOST_UBLAS_CHECK (filled_ <= non_zeros, internal_logic ()); } // Proxy support @@ -4135,7 +4134,22 @@ namespace boost { namespace numeric { namespace ublas { index_triple_array ita (filled_, index1_data_, index2_data_, value_data_); std::sort (ita.begin (), ita.end ()); - filled_ = std::unique (ita.begin (), ita.end ()) - ita.begin(); + // ISSUE: unusual semantics - sum values of duplicates + size_type filled = 1; + for (size_type i = 1; i < filled_; ++ i) { + if (index1_data_ [filled - 1] != index1_data_ [i] || + index2_data_ [filled - 1] != index2_data_ [i]) { + ++ filled; + if (filled - 1 != i) { + index1_data_ [filled - 1] = index1_data_ [i]; + index2_data_ [filled - 1] = index2_data_ [i]; + value_data_ [filled - 1] = value_data_ [i]; + } + } else { + value_data_ [filled - 1] += value_data_ [i]; + } + } + filled_ = filled; sorted_ = true; } } @@ -4143,17 +4157,18 @@ namespace boost { namespace numeric { namespace ublas { // Element insertion and erasure BOOST_UBLAS_INLINE void push_back (size_type i, size_type j, const_reference t) { - if (filled_ >= non_zeros_) - reserve (2 * non_zeros_, true); size_type element1 = functor_type::element1 (i, size1_, j, size2_); size_type element2 = functor_type::element2 (i, size1_, j, size2_); if (filled_ == 0 || index1_data () [filled_ - 1] < k_based (element1) || (index1_data () [filled_ - 1] == k_based (element1) && index2_data () [filled_ - 1] < k_based (element2))) { + if (filled_ >= non_zeros_) + reserve (2 * filled_, true); + BOOST_UBLAS_CHECK (filled_ < non_zeros_, internal_logic ()); + index1_data () [filled_] = k_based (element1); + index2_data () [filled_] = k_based (element2); + value_data () [filled_] = t; ++ filled_; - index1_data () [filled_ - 1] = k_based (element1); - index2_data () [filled_ - 1] = k_based (element2); - value_data () [filled_ - 1] = t; return; } external_logic ().raise (); @@ -4161,13 +4176,14 @@ namespace boost { namespace numeric { namespace ublas { BOOST_UBLAS_INLINE void insert (size_type i, size_type j, const_reference t) { if (filled_ >= non_zeros_) - reserve (2 * non_zeros_, true); + reserve (2 * filled_, true); + BOOST_UBLAS_CHECK (filled_ < non_zeros_, internal_logic ()); size_type element1 = functor_type::element1 (i, size1_, j, size2_); size_type element2 = functor_type::element2 (i, size1_, j, size2_); + index1_data () [filled_] = k_based (element1); + index2_data () [filled_] = k_based (element2); + value_data () [filled_] = t; ++ filled_; - index1_data () [filled_ - 1] = k_based (element1); - index2_data () [filled_ - 1] = k_based (element2); - value_data () [filled_ - 1] = t; sorted_ = false; } BOOST_UBLAS_INLINE diff --git a/include/boost/numeric/ublas/vector_sparse.hpp b/include/boost/numeric/ublas/vector_sparse.hpp index c265b15f..f62e4b4a 100644 --- a/include/boost/numeric/ublas/vector_sparse.hpp +++ b/include/boost/numeric/ublas/vector_sparse.hpp @@ -720,12 +720,12 @@ namespace boost { namespace numeric { namespace ublas { BOOST_UBLAS_INLINE compressed_vector (): vector_expression (), - size_ (0), non_zeros_ (max_nz (0)), filled_ (0), + size_ (0), non_zeros_ (restrict_nz (0)), filled_ (0), index_data_ (non_zeros_), value_data_ (non_zeros_) {} explicit BOOST_UBLAS_INLINE compressed_vector (size_type size, size_type non_zeros = 0): vector_expression (), - size_ (size), non_zeros_ (max_nz (non_zeros)), filled_ (0), + size_ (size), non_zeros_ (restrict_nz (non_zeros)), filled_ (0), index_data_ (non_zeros_), value_data_ (non_zeros_) { } BOOST_UBLAS_INLINE @@ -737,7 +737,7 @@ namespace boost { namespace numeric { namespace ublas { BOOST_UBLAS_INLINE compressed_vector (const vector_expression &ae, size_type non_zeros = 0): vector_expression (), - size_ (ae ().size ()), non_zeros_ (max_nz (non_zeros)), filled_ (0), + size_ (ae ().size ()), non_zeros_ (restrict_nz (non_zeros)), filled_ (0), index_data_ (non_zeros_), value_data_ (non_zeros_) { vector_assign (scalar_assign (), *this, ae); } @@ -779,7 +779,7 @@ namespace boost { namespace numeric { namespace ublas { // Resizing private: BOOST_UBLAS_INLINE - size_type max_nz (size_type non_zeros) const { + size_type restrict_nz (size_type non_zeros) const { non_zeros = (std::max) (non_zeros, size_type (1)); non_zeros = (std::min) (non_zeros, size_); return non_zeros; @@ -790,7 +790,7 @@ namespace boost { namespace numeric { namespace ublas { // FIXME preserve unimplemented BOOST_UBLAS_CHECK (!preserve, internal_logic ()); size_ = size; - non_zeros_ = max_nz (non_zeros_); + non_zeros_ = restrict_nz (non_zeros_); index_data (). resize (non_zeros_); value_data (). resize (non_zeros_); filled_ = 0; @@ -799,7 +799,7 @@ namespace boost { namespace numeric { namespace ublas { // Reserving BOOST_UBLAS_INLINE void reserve (size_type non_zeros, bool preserve = true) { - non_zeros_ = max_nz (non_zeros); + non_zeros_ = restrict_nz (non_zeros); if (preserve) { index_data (). resize (non_zeros_, size_type ()); value_data (). resize (non_zeros_, value_type ()); @@ -1280,12 +1280,12 @@ namespace boost { namespace numeric { namespace ublas { BOOST_UBLAS_INLINE coordinate_vector (): vector_expression (), - size_ (0), non_zeros_ (max_nz (0)), filled_ (0), + size_ (0), non_zeros_ (restrict_nz (0)), filled_ (0), sorted_ (true), index_data_ (non_zeros_), value_data_ (non_zeros_) {} explicit BOOST_UBLAS_INLINE coordinate_vector (size_type size, size_type non_zeros = 0): vector_expression (), - size_ (size), non_zeros_ (max_nz (non_zeros)), filled_ (0), + size_ (size), non_zeros_ (restrict_nz (non_zeros)), filled_ (0), sorted_ (true), index_data_ (non_zeros_), value_data_ (non_zeros_) { } BOOST_UBLAS_INLINE @@ -1297,7 +1297,7 @@ namespace boost { namespace numeric { namespace ublas { BOOST_UBLAS_INLINE coordinate_vector (const vector_expression &ae, size_type non_zeros = 0): vector_expression (), - size_ (ae ().size ()), non_zeros_ (max_nz (non_zeros)), filled_ (0), + size_ (ae ().size ()), non_zeros_ (restrict_nz (non_zeros)), filled_ (0), sorted_ (true), index_data_ (non_zeros_), value_data_ (non_zeros_) { vector_assign (scalar_assign (), *this, ae); } @@ -1339,19 +1339,18 @@ namespace boost { namespace numeric { namespace ublas { // Resizing private: BOOST_UBLAS_INLINE - size_type max_nz (size_type non_zeros) const { + size_type restrict_nz (size_type non_zeros) const { + // minimum non_zeros non_zeros = (std::max) (non_zeros, size_type (1)); - non_zeros = (std::min) (non_zeros, size_); + // ISSUE no maximum as coordinate may contain inserted duplicates return non_zeros; } public: BOOST_UBLAS_INLINE void resize (size_type size, bool preserve = true) { - // FIXME preserve unimplemented - BOOST_UBLAS_CHECK (!preserve, internal_logic ()); if (preserve) sort (); // remove duplicate elements. - non_zeros_ = max_nz (non_zeros_); + non_zeros_ = restrict_nz (non_zeros_); if (preserve) { index_data (). resize (non_zeros_, size_type ()); value_data (). resize (non_zeros_, value_type ()); @@ -1363,13 +1362,14 @@ namespace boost { namespace numeric { namespace ublas { filled_ = 0; } size_ = size; + BOOST_UBLAS_CHECK (filled_ <= non_zeros, internal_logic ()); } // Reserving BOOST_UBLAS_INLINE void reserve (size_type non_zeros, bool preserve = true) { if (preserve) sort (); // remove duplicate elements. - non_zeros_ = max_nz (non_zeros); + non_zeros_ = restrict_nz (non_zeros); if (preserve) { index_data (). resize (non_zeros_, size_type ()); value_data (). resize (non_zeros_, value_type ()); @@ -1380,6 +1380,7 @@ namespace boost { namespace numeric { namespace ublas { value_data (). resize (non_zeros_); filled_ = 0; } + BOOST_UBLAS_CHECK (filled_ <= non_zeros, internal_logic ()); } // Proxy support @@ -1527,7 +1528,7 @@ namespace boost { namespace numeric { namespace ublas { index_pair_array ipa (filled_, index_data_, value_data_); std::sort (ipa.begin (), ipa.end ()); - // FIX: check for duplicates + // ISSUE: unusual semantics - sum values of duplicates size_type filled = 1; for (size_type i = 1; i < filled_; ++ i) { if (index_data_ [filled - 1] != index_data_ [i]) { @@ -1548,12 +1549,13 @@ namespace boost { namespace numeric { namespace ublas { // Element insertion and erasure BOOST_UBLAS_INLINE void push_back (size_type i, const_reference t) { - if (filled_ >= non_zeros_) - reserve (2 * non_zeros_, true); if (filled_ == 0 || index_data () [filled_ - 1] < k_based (i)) { + if (filled_ >= non_zeros_) + reserve (2 * filled_, true); + BOOST_UBLAS_CHECK (filled_ < non_zeros_, internal_logic ()); + index_data () [filled_] = k_based (i); + value_data () [filled_] = t; ++ filled_; - index_data () [filled_ - 1] = k_based (i); - value_data () [filled_ - 1] = t; return; } external_logic ().raise (); @@ -1561,10 +1563,11 @@ namespace boost { namespace numeric { namespace ublas { BOOST_UBLAS_INLINE void insert (size_type i, const_reference t) { if (filled_ >= non_zeros_) - reserve (2 * non_zeros_, true); + reserve (2 * filled_, true); + BOOST_UBLAS_CHECK (filled_ < non_zeros_, internal_logic ()); + index_data () [filled_] = k_based (i); + value_data () [filled_] = t; ++ filled_; - index_data () [filled_ - 1] = k_based (i); - value_data () [filled_ - 1] = t; sorted_ = false; } BOOST_UBLAS_INLINE