Index: boost/numeric/ublas/vector_sparse.hpp =================================================================== --- boost/numeric/ublas/vector_sparse.hpp (revision 80623) +++ boost/numeric/ublas/vector_sparse.hpp (working copy) @@ -262,15 +262,15 @@ /** \brief Index map based sparse vector * - * A sparse vector of values of type T of variable size. The sparse storage type A can be + * A sparse vector of values of type T of variable size. The sparse storage type A can be * \c std::map or \c map_array. This means that only non-zero elements * are effectively stored. * - * For a \f$n\f$-dimensional sparse vector, and 0 <= i < n the non-zero elements \f$v_i\f$ - * are mapped to consecutive elements of the associative container, i.e. for elements + * For a \f$n\f$-dimensional sparse vector, and 0 <= i < n the non-zero elements \f$v_i\f$ + * are mapped to consecutive elements of the associative container, i.e. for elements * \f$k = v_{i_1}\f$ and \f$k + 1 = v_{i_2}\f$ of the container, holds \f$i_1 < i_2\f$. * - * Supported parameters for the adapted array are \c map_array and + * Supported parameters for the adapted array are \c map_array and * \c map_std. The latter is equivalent to \c std::map. * * \tparam T the type of object stored in the vector (like double, float, complex, etc...) @@ -771,18 +771,18 @@ // Thanks to Kresimir Fresl for extending this to cover different index bases. - + /** \brief Compressed array based sparse vector * - * a sparse vector of values of type T of variable size. The non zero values are stored as - * two seperate arrays: an index array and a value array. The index array is always sorted + * a sparse vector of values of type T of variable size. The non zero values are stored as + * two seperate arrays: an index array and a value array. The index array is always sorted * and there is at most one entry for each index. Inserting an element can be time consuming. * If the vector contains a few zero entries, then it is better to have a normal vector. * If the vector has a very high dimension with a few non-zero values, then this vector is * very memory efficient (at the cost of a few more computations). * - * For a \f$n\f$-dimensional compressed vector and \f$0 \leq i < n\f$ the non-zero elements - * \f$v_i\f$ are mapped to consecutive elements of the index and value container, i.e. for + * For a \f$n\f$-dimensional compressed vector and \f$0 \leq i < n\f$ the non-zero elements + * \f$v_i\f$ are mapped to consecutive elements of the index and value container, i.e. for * elements \f$k = v_{i_1}\f$ and \f$k + 1 = v_{i_2}\f$ of these containers holds \f$i_1 < i_2\f$. * * Supported parameters for the adapted array (indices and values) are \c unbounded_array<> , @@ -1414,15 +1414,15 @@ /** \brief Coordimate array based sparse vector * - * a sparse vector of values of type \c T of variable size. The non zero values are stored - * as two seperate arrays: an index array and a value array. The arrays may be out of order - * with multiple entries for each vector element. If there are multiple values for the same + * a sparse vector of values of type \c T of variable size. The non zero values are stored + * as two seperate arrays: an index array and a value array. The arrays may be out of order + * with multiple entries for each vector element. If there are multiple values for the same * index the sum of these values is the real value. It is way more efficient for inserting values - * than a \c compressed_vector but less memory efficient. Also linearly parsing a vector can + * than a \c compressed_vector but less memory efficient. Also linearly parsing a vector can * be longer in specific cases than a \c compressed_vector. * - * For a n-dimensional sorted coordinate vector and \f$ 0 \leq i < n\f$ the non-zero elements - * \f$v_i\f$ are mapped to consecutive elements of the index and value container, i.e. for + * For a n-dimensional sorted coordinate vector and \f$ 0 \leq i < n\f$ the non-zero elements + * \f$v_i\f$ are mapped to consecutive elements of the index and value container, i.e. for * elements \f$k = v_{i_1}\f$ and \f$k + 1 = v_{i_2}\f$ of these containers holds \f$i_1 < i_2\f$. * * Supported parameters for the adapted array (indices and values) are \c unbounded_array<> , @@ -1801,6 +1801,48 @@ v1.swap (v2); } + // replacement if STL lower bound algorithm for use of inplace_merge + size_type lower_bound (size_type beg, size_type end, size_type target) const { + while (end > beg) { + size_type mid = (beg + end) / 2; + if (index_data_[mid] < index_data_[target]) { + beg = mid + 1; + } else { + end = mid; + } + } + return beg; + } + + // specialized replacement of STL inplace_merge to avoid compilation + // problems with respect to the array_triple iterator + void inplace_merge (size_type beg, size_type mid, size_type end) const { + size_type len_lef = mid - beg; + size_type len_rig = end - mid; + + if (len_lef == 1 && len_rig == 1) { + if (index_data_[mid] < index_data_[beg]) { + std::swap(index_data_[beg], index_data_[mid]); + std::swap(value_data_[beg], value_data_[mid]); + } + } else if (len_lef > 0 && len_rig > 0) { + size_type lef_mid, rig_mid; + if (len_lef >= len_rig) { + lef_mid = (beg + mid) / 2; + rig_mid = lower_bound(mid, end, lef_mid); + } else { + rig_mid = (mid + end) / 2; + lef_mid = lower_bound(beg, mid, rig_mid); + } + std::rotate(&index_data_[0] + lef_mid, &index_data_[0] + mid, &index_data_[0] + rig_mid); + std::rotate(&value_data_[0] + lef_mid, &value_data_[0] + mid, &value_data_[0] + rig_mid); + + size_type new_mid = lef_mid + rig_mid - mid; + inplace_merge(beg, lef_mid, new_mid); + inplace_merge(new_mid, rig_mid, end); + } + } + // Sorting and summation of duplicates BOOST_UBLAS_INLINE void sort () const { @@ -1811,11 +1853,11 @@ const typename array_pair::iterator iunsorted = ipa.begin () + sorted_filled_; // sort new elements and merge std::sort (iunsorted, ipa.end ()); - std::inplace_merge (ipa.begin (), iunsorted, ipa.end ()); + inplace_merge(0, sorted_filled_, filled_); #else const typename array_pair::iterator iunsorted = ipa.begin (); std::sort (iunsorted, ipa.end ()); -#endif +#endif // sum duplicates with += and remove size_type filled = 0; Index: boost/numeric/ublas/matrix_sparse.hpp =================================================================== --- boost/numeric/ublas/matrix_sparse.hpp (revision 80623) +++ boost/numeric/ublas/matrix_sparse.hpp (working copy) @@ -54,7 +54,7 @@ else *p = s; } - + public: // Construction and destruction BOOST_UBLAS_INLINE @@ -250,14 +250,14 @@ * * This class represents a matrix by using a \c key to value mapping. The default type is * \code template > class mapped_matrix; \endcode - * So, by default a STL map container is used to associate keys and values. The key is computed depending on + * So, by default a STL map container is used to associate keys and values. The key is computed depending on * the layout type \c L as \code key = layout_type::element(i, size1_, j, size2_); \endcode * which means \code key = (i*size2+j) \endcode for a row major matrix. - * Limitations: The matrix size must not exceed \f$(size1*size2) < \f$ \code std::limits \endcode. + * Limitations: The matrix size must not exceed \f$(size1*size2) < \f$ \code std::limits \endcode. * The \ref find1() and \ref find2() operations have a complexity of at least \f$\mathcal{O}(log(nnz))\f$, depending * on the efficiency of \c std::lower_bound on the key set of the map. - * Orientation and storage can also be specified, otherwise a row major orientation is used. - * It is \b not required by the storage to initialize elements of the matrix. By default, the orientation is \c row_major. + * Orientation and storage can also be specified, otherwise a row major orientation is used. + * It is \b not required by the storage to initialize elements of the matrix. By default, the orientation is \c row_major. * * \sa fwd.hpp, storage_sparse.hpp * @@ -427,7 +427,7 @@ return; data ().erase (it); } - + // Zeroing BOOST_UBLAS_INLINE void clear () { @@ -557,7 +557,7 @@ typedef reverse_iterator_base2 reverse_iterator2; // Element lookup - // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. + // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const { const_subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_))); const_subiterator_type it_end (data ().end ()); @@ -586,7 +586,7 @@ } return const_iterator1 (*this, rank, i, j, it); } - // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. + // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) { subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_))); subiterator_type it_end (data ().end ()); @@ -615,7 +615,7 @@ } return iterator1 (*this, rank, i, j, it); } - // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. + // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const { const_subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_))); const_subiterator_type it_end (data ().end ()); @@ -644,7 +644,7 @@ } return const_iterator2 (*this, rank, i, j, it); } - // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. + // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) { subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_))); subiterator_type it_end (data ().end ()); @@ -1509,7 +1509,7 @@ return; (*itv).second.erase (it); } - + // Zeroing BOOST_UBLAS_INLINE void clear () { @@ -1632,7 +1632,7 @@ subiterator_type it ((*itv).second.find (element2)); BOOST_UBLAS_CHECK (it != (*itv).second.end (), bad_index ()); BOOST_UBLAS_CHECK ((*it).first == element2, internal_logic ()); // broken map - + return it->second; } @@ -1647,7 +1647,7 @@ typedef reverse_iterator_base2 reverse_iterator2; // Element lookup - // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. + // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const { BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ()); for (;;) { @@ -1662,8 +1662,8 @@ // advance to the first available major index size_type M = itv->first; size_type m; - if (it != it_end) { - m = it->first; + if (it != it_end) { + m = it->first; } else { m = layout_type::size_m(size1_, size2_); } @@ -1696,7 +1696,7 @@ } } } - // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. + // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) { BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ()); for (;;) { @@ -1711,8 +1711,8 @@ // advance to the first available major index size_type M = itv->first; size_type m; - if (it != it_end) { - m = it->first; + if (it != it_end) { + m = it->first; } else { m = layout_type::size_m(size1_, size2_); } @@ -1745,7 +1745,7 @@ } } } - // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. + // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const { BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ()); for (;;) { @@ -1760,8 +1760,8 @@ // advance to the first available major index size_type M = itv->first; size_type m; - if (it != it_end) { - m = it->first; + if (it != it_end) { + m = it->first; } else { m = layout_type::size_m(size1_, size2_); } @@ -1794,7 +1794,7 @@ } } } - // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. + // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) { BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ()); for (;;) { @@ -1809,8 +1809,8 @@ // advance to the first available major index size_type M = itv->first; size_type m; - if (it != it_end) { - m = it->first; + if (it != it_end) { + m = it->first; } else { m = layout_type::size_m(size1_, size2_); } @@ -2682,7 +2682,7 @@ index1_data_ (m.index1_data_), index2_data_ (m.index2_data_), value_data_ (m.value_data_) { storage_invariants (); } - + BOOST_UBLAS_INLINE compressed_matrix (const coordinate_matrix &m): matrix_container (), @@ -2737,7 +2737,7 @@ size_type nnz () const { return filled2_; } - + // Storage accessors BOOST_UBLAS_INLINE static size_type index_base () { @@ -2942,7 +2942,7 @@ } storage_invariants (); } - + // Zeroing BOOST_UBLAS_INLINE void clear () { @@ -3121,7 +3121,7 @@ typedef reverse_iterator_base2 reverse_iterator2; // Element lookup - // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. + // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const { for (;;) { array_size_type address1 (layout_type::index_M (i, j)); @@ -3161,7 +3161,7 @@ } } } - // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. + // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) { for (;;) { array_size_type address1 (layout_type::index_M (i, j)); @@ -3201,7 +3201,7 @@ } } } - // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. + // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const { for (;;) { array_size_type address1 (layout_type::index_M (i, j)); @@ -3241,7 +3241,7 @@ } } } - // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. + // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) { for (;;) { array_size_type address1 (layout_type::index_M (i, j)); @@ -3968,7 +3968,7 @@ BOOST_UBLAS_CHECK (filled2_ <= capacity_, internal_logic ()); BOOST_UBLAS_CHECK (index1_data_ [filled1_ - 1] == k_based (filled2_), internal_logic ()); } - + size_type size1_; size_type size2_; array_size_type capacity_; @@ -4208,7 +4208,7 @@ const_pointer p = find_element (i, j); if (p) return *p; - else + else return zero_; } BOOST_UBLAS_INLINE @@ -4267,7 +4267,7 @@ } storage_invariants (); } - + // Zeroing BOOST_UBLAS_INLINE void clear () { @@ -4391,6 +4391,54 @@ m1.swap (m2); } + // replacement if STL lower bound algorithm for use of inplace_merge + array_size_type lower_bound (array_size_type beg, array_size_type end, array_size_type target) const { + while (end > beg) { + array_size_type mid = (beg + end) / 2; + if (((index1_data_[mid] < index1_data_[target]) || + ((index1_data_[mid] == index1_data_[target]) && + (index2_data_[mid] < index2_data_[target])))) { + beg = mid + 1; + } else { + end = mid; + } + } + return beg; + } + + // specialized replacement of STL inplace_merge to avoid compilation + // problems with respect to the array_triple iterator + void inplace_merge (array_size_type beg, array_size_type mid, array_size_type end) const { + array_size_type len_lef = mid - beg; + array_size_type len_rig = end - mid; + + if (len_lef == 1 && len_rig == 1) { + if ((index1_data_[mid] < index1_data_[beg]) || + ((index1_data_[mid] == index1_data_[beg]) && (index2_data_[mid] < index2_data_[beg]))) + { + std::swap(index1_data_[beg], index1_data_[mid]); + std::swap(index2_data_[beg], index2_data_[mid]); + std::swap(value_data_[beg], value_data_[mid]); + } + } else if (len_lef > 0 && len_rig > 0) { + array_size_type lef_mid, rig_mid; + if (len_lef >= len_rig) { + lef_mid = (beg + mid) / 2; + rig_mid = lower_bound(mid, end, lef_mid); + } else { + rig_mid = (mid + end) / 2; + lef_mid = lower_bound(beg, mid, rig_mid); + } + std::rotate(&index1_data_[0] + lef_mid, &index1_data_[0] + mid, &index1_data_[0] + rig_mid); + std::rotate(&index2_data_[0] + lef_mid, &index2_data_[0] + mid, &index2_data_[0] + rig_mid); + std::rotate(&value_data_[0] + lef_mid, &value_data_[0] + mid, &value_data_[0] + rig_mid); + + array_size_type new_mid = lef_mid + rig_mid - mid; + inplace_merge(beg, lef_mid, new_mid); + inplace_merge(new_mid, rig_mid, end); + } + } + // Sorting and summation of duplicates BOOST_UBLAS_INLINE void sort () const { @@ -4401,11 +4449,11 @@ const typename array_triple::iterator iunsorted = ita.begin () + sorted_filled_; // sort new elements and merge std::sort (iunsorted, ita.end ()); - std::inplace_merge (ita.begin (), iunsorted, ita.end ()); + inplace_merge(0, sorted_filled_, filled_); #else const typename array_triple::iterator iunsorted = ita.begin (); std::sort (iunsorted, ita.end ()); -#endif +#endif // sum duplicates with += and remove array_size_type filled = 0; for (array_size_type i = 1; i < filled_; ++ i) { @@ -4434,7 +4482,7 @@ size_type element1 = layout_type::index_M (i, j); size_type element2 = layout_type::index_m (i, j); // must maintain sort order - BOOST_UBLAS_CHECK (sorted_ && + BOOST_UBLAS_CHECK (sorted_ && (filled_ == 0 || index1_data_ [filled_ - 1] < k_based (element1) || (index1_data_ [filled_ - 1] == k_based (element1) && index2_data_ [filled_ - 1] < k_based (element2))) @@ -4485,7 +4533,7 @@ typedef reverse_iterator_base2 reverse_iterator2; // Element lookup - // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. + // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const { sort (); for (;;) { @@ -4526,7 +4574,7 @@ } } } - // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. + // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) { sort (); for (;;) { @@ -4567,7 +4615,7 @@ } } } - // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. + // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const { sort (); for (;;) { @@ -4608,7 +4656,7 @@ } } } - // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. + // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) { sort (); for (;;) { Index: libs/numeric/ublas/test/test_coordinate_vector_inplace_merge.cpp =================================================================== --- libs/numeric/ublas/test/test_coordinate_vector_inplace_merge.cpp (revision 0) +++ libs/numeric/ublas/test/test_coordinate_vector_inplace_merge.cpp (revision 0) @@ -0,0 +1,120 @@ +// Copyright (c) 2011 David Bellot +// +// Distributed under the Boost Software License, Version 1.0. (See +// accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_UBLAS_NO_ELEMENT_PROXIES +# define BOOST_UBLAS_NO_ELEMENT_PROXIES +#endif + +#include +#include +#include +#include +#include + +#include "libs/numeric/ublas/test/utils.hpp" + +const double TOL = 1e-15; + + +template +typename AE::value_type mean_square(const boost::numeric::ublas::vector_expression &me) { + typename AE::value_type s(0); + typename AE::size_type i; + for (i=0; i!= me().size(); i++) { + s += boost::numeric::ublas::scalar_traits::type_abs(me()(i)); + } + return s/me().size(); +} + +template +bool check_sortedness(const boost::numeric::ublas::coordinate_vector& vector) { + bool result = true; + typedef boost::numeric::ublas::coordinate_vector vector_type; + typename vector_type::index_array_type idx = vector.index_data(); + typename vector_type::size_type size = vector.filled(); + + for (typename vector_type::size_type i = 0; i + 1 < size && result; ++ i) { + result &= (idx[i] < idx[i + 1]); + } + return result; +} + +void print_entries(size_t size, + const std::vector& entries) +{ + std::cerr << "Error entries - Size:" << size << ". Entries: "; + for (size_t i = 0; i < entries.size(); ++ i) { + std::cerr << entries[i] << "; "; + } + std::cerr << "\n"; +} + +BOOST_UBLAS_TEST_DEF( test_coordinate_vector_inplace_merge_random ) +{ + const size_t max_repeats = 100; + const size_t max_size = 100; + const size_t dim_var = 10; + const size_t nr_entries = 10; + + for (size_t repeats = 1; repeats < max_repeats; ++repeats ) { + for (size_t size = 1; size < max_size; size += 5) { + size_t size_vec = size + rand() % dim_var; + + boost::numeric::ublas::coordinate_vector vector_coord(size_vec); + boost::numeric::ublas::vector vector_dense(size_vec, 0); + + vector_coord.sort(); + + std::vector entries; + for (int entry = 0; entry < nr_entries; ++ entry) { + int x = rand() % size_vec; + entries.push_back(x); + vector_coord.append_element(x, 1); + vector_dense(x) += 1; + } + vector_coord.sort(); + + { + bool sorted = check_sortedness(vector_coord); + bool identical = mean_square(vector_coord - vector_dense) < TOL; + if (!(sorted && identical)) { + print_entries(size_vec, entries); + } + BOOST_UBLAS_TEST_CHECK( check_sortedness(vector_coord) ); + BOOST_UBLAS_TEST_CHECK( mean_square(vector_coord - vector_dense) < TOL); + } + + for (int entry = 0; entry < nr_entries; ++ entry) { + int x = rand() % size_vec; + entries.push_back(x); + vector_coord(x) += 1; + vector_dense(x) += 1; + vector_coord.sort(); + } + + { + bool sorted = check_sortedness(vector_coord); + bool identical = mean_square(vector_coord - vector_dense) < TOL; + if (!(sorted && identical)) { + print_entries(size_vec, entries); + } + BOOST_UBLAS_TEST_CHECK( sorted ); + BOOST_UBLAS_TEST_CHECK( identical ); + } + } + } +} + +int main() +{ + BOOST_UBLAS_TEST_BEGIN(); + + BOOST_UBLAS_TEST_DO( test_coordinate_vector_inplace_merge_random ); + + BOOST_UBLAS_TEST_END(); + + return EXIT_SUCCESS;; +} Index: libs/numeric/ublas/test/test_coordinate_matrix_inplace_merge.cpp =================================================================== --- libs/numeric/ublas/test/test_coordinate_matrix_inplace_merge.cpp (revision 0) +++ libs/numeric/ublas/test/test_coordinate_matrix_inplace_merge.cpp (revision 0) @@ -0,0 +1,132 @@ +// Copyright (c) 2011 David Bellot +// +// Distributed under the Boost Software License, Version 1.0. (See +// accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_UBLAS_NO_ELEMENT_PROXIES +# define BOOST_UBLAS_NO_ELEMENT_PROXIES +#endif + +#include +#include +#include +#include +#include + +#include "libs/numeric/ublas/test/utils.hpp" + +using std::cout; +using std::endl; + +const double TOL = 1e-15; + +template +typename AE::value_type mean_square(const boost::numeric::ublas::matrix_expression &me) { + typename AE::value_type s(0); + typename AE::size_type i, j; + for (i=0; i!= me().size1(); i++) { + for (j=0; j!= me().size2(); j++) { + s += boost::numeric::ublas::scalar_traits::type_abs(me()(i,j)); + } + } + return s/me().size1()*me().size2(); +} + +template +bool check_sortedness(const boost::numeric::ublas::coordinate_matrix& matrix) { + bool result = true; + typedef boost::numeric::ublas::coordinate_matrix matrix_type; + typename matrix_type::index_array_type i1 = matrix.index1_data(); + typename matrix_type::index_array_type i2 = matrix.index2_data(); + typename matrix_type::array_size_type size = matrix.filled(); + + for (typename matrix_type::array_size_type i = 0; i + 1 < size && result; ++ i) { + result &= ( (i1[i] < i1[i + 1]) || + ((i1[i] == i1[i]) && + (i2[i] < i2[i + 1])) ); + + } + return result; +} + +void print_entries(size_t size_x, size_t size_y, + const std::vector >& entries) +{ + std::cerr << "Error - Size:" << size_x << " x " << size_y << ". Entries: "; + for (size_t i = 0; i < entries.size(); ++ i) { + std::cerr << entries[i].first << ", " << entries[i].second << "; "; + } + std::cerr << "\n"; +} + + +BOOST_UBLAS_TEST_DEF( test_coordinate_matrix_inplace_merge_random ) +{ + const size_t max_repeats = 100; + const size_t max_size = 100; + const size_t dim_var = 10; + const size_t nr_entries = 10; + + for (size_t repeats = 1; repeats < max_repeats; ++repeats ) { + for (size_t size = 1; size < max_size; size += 5) { + size_t size_x = size + rand() % dim_var; + size_t size_y = size + rand() % dim_var; + + boost::numeric::ublas::coordinate_matrix matrix_coord(size_x, size_y); + boost::numeric::ublas::matrix matrix_dense(size_x, size_y, 0); + + matrix_coord.sort(); + + std::vector > entries; + for (int entry = 0; entry < nr_entries; ++ entry) { + int x = rand() % size_x; + int y = rand() % size_y; + entries.push_back(std::make_pair(x, y)); + matrix_coord.append_element(x, y, 1); + matrix_dense(x, y) += 1; + } + matrix_coord.sort(); + + { + bool sorted = check_sortedness(matrix_coord); + bool identical = mean_square(matrix_coord - matrix_dense) < TOL; + if (!(sorted && identical)) { + print_entries(size_x, size_y, entries); + } + BOOST_UBLAS_TEST_CHECK( check_sortedness(matrix_coord) ); + BOOST_UBLAS_TEST_CHECK( mean_square(matrix_coord - matrix_dense) < TOL); + } + + for (int entry = 0; entry < nr_entries; ++ entry) { + int x = rand() % size_x; + int y = rand() % size_y; + entries.push_back(std::make_pair(x, y)); + matrix_coord(x, y) += 1; + matrix_dense(x, y) += 1; + matrix_coord.sort(); + } + + { + bool sorted = check_sortedness(matrix_coord); + bool identical = mean_square(matrix_coord - matrix_dense) < TOL; + if (!(sorted && identical)) { + print_entries(size_x, size_y, entries); + } + BOOST_UBLAS_TEST_CHECK( sorted ); + BOOST_UBLAS_TEST_CHECK( identical ); + } + } + } +} + +int main() +{ + BOOST_UBLAS_TEST_BEGIN(); + + BOOST_UBLAS_TEST_DO( test_coordinate_matrix_inplace_merge_random ); + + BOOST_UBLAS_TEST_END(); + + return EXIT_SUCCESS;; +}