Ticket #7363: custom_inplace_merge_fixed.diff
| File custom_inplace_merge_fixed.diff, 33.4 KB (added by , 10 years ago) |
|---|
-
boost/numeric/ublas/vector_sparse.hpp
262 262 263 263 /** \brief Index map based sparse vector 264 264 * 265 * A sparse vector of values of type T of variable size. The sparse storage type A can be 265 * A sparse vector of values of type T of variable size. The sparse storage type A can be 266 266 * \c std::map<size_t, T> or \c map_array<size_t, T>. This means that only non-zero elements 267 267 * are effectively stored. 268 268 * 269 * For a \f$n\f$-dimensional sparse vector, and 0 <= i < n the non-zero elements \f$v_i\f$ 270 * are mapped to consecutive elements of the associative container, i.e. for elements 269 * For a \f$n\f$-dimensional sparse vector, and 0 <= i < n the non-zero elements \f$v_i\f$ 270 * are mapped to consecutive elements of the associative container, i.e. for elements 271 271 * \f$k = v_{i_1}\f$ and \f$k + 1 = v_{i_2}\f$ of the container, holds \f$i_1 < i_2\f$. 272 272 * 273 * Supported parameters for the adapted array are \c map_array<std::size_t, T> and 273 * Supported parameters for the adapted array are \c map_array<std::size_t, T> and 274 274 * \c map_std<std::size_t, T>. The latter is equivalent to \c std::map<std::size_t, T>. 275 275 * 276 276 * \tparam T the type of object stored in the vector (like double, float, complex, etc...) … … 771 771 772 772 773 773 // Thanks to Kresimir Fresl for extending this to cover different index bases. 774 774 775 775 /** \brief Compressed array based sparse vector 776 776 * 777 * a sparse vector of values of type T of variable size. The non zero values are stored as 778 * two seperate arrays: an index array and a value array. The index array is always sorted 777 * a sparse vector of values of type T of variable size. The non zero values are stored as 778 * two seperate arrays: an index array and a value array. The index array is always sorted 779 779 * and there is at most one entry for each index. Inserting an element can be time consuming. 780 780 * If the vector contains a few zero entries, then it is better to have a normal vector. 781 781 * If the vector has a very high dimension with a few non-zero values, then this vector is 782 782 * very memory efficient (at the cost of a few more computations). 783 783 * 784 * For a \f$n\f$-dimensional compressed vector and \f$0 \leq i < n\f$ the non-zero elements 785 * \f$v_i\f$ are mapped to consecutive elements of the index and value container, i.e. for 784 * For a \f$n\f$-dimensional compressed vector and \f$0 \leq i < n\f$ the non-zero elements 785 * \f$v_i\f$ are mapped to consecutive elements of the index and value container, i.e. for 786 786 * 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$. 787 787 * 788 788 * Supported parameters for the adapted array (indices and values) are \c unbounded_array<> , … … 1414 1414 1415 1415 /** \brief Coordimate array based sparse vector 1416 1416 * 1417 * a sparse vector of values of type \c T of variable size. The non zero values are stored 1418 * as two seperate arrays: an index array and a value array. The arrays may be out of order 1419 * with multiple entries for each vector element. If there are multiple values for the same 1417 * a sparse vector of values of type \c T of variable size. The non zero values are stored 1418 * as two seperate arrays: an index array and a value array. The arrays may be out of order 1419 * with multiple entries for each vector element. If there are multiple values for the same 1420 1420 * index the sum of these values is the real value. It is way more efficient for inserting values 1421 * than a \c compressed_vector but less memory efficient. Also linearly parsing a vector can 1421 * than a \c compressed_vector but less memory efficient. Also linearly parsing a vector can 1422 1422 * be longer in specific cases than a \c compressed_vector. 1423 1423 * 1424 * For a n-dimensional sorted coordinate vector and \f$ 0 \leq i < n\f$ the non-zero elements 1425 * \f$v_i\f$ are mapped to consecutive elements of the index and value container, i.e. for 1424 * For a n-dimensional sorted coordinate vector and \f$ 0 \leq i < n\f$ the non-zero elements 1425 * \f$v_i\f$ are mapped to consecutive elements of the index and value container, i.e. for 1426 1426 * 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$. 1427 1427 * 1428 1428 * Supported parameters for the adapted array (indices and values) are \c unbounded_array<> , … … 1801 1801 v1.swap (v2); 1802 1802 } 1803 1803 1804 // replacement if STL lower bound algorithm for use of inplace_merge 1805 size_type lower_bound (size_type beg, size_type end, size_type target) const { 1806 while (end > beg) { 1807 size_type mid = (beg + end) / 2; 1808 if (index_data_[mid] < index_data_[target]) { 1809 beg = mid + 1; 1810 } else { 1811 end = mid; 1812 } 1813 } 1814 return beg; 1815 } 1816 1817 // specialized replacement of STL inplace_merge to avoid compilation 1818 // problems with respect to the array_triple iterator 1819 void inplace_merge (size_type beg, size_type mid, size_type end) const { 1820 size_type len_lef = mid - beg; 1821 size_type len_rig = end - mid; 1822 1823 if (len_lef == 1 && len_rig == 1) { 1824 if (index_data_[mid] < index_data_[beg]) { 1825 std::swap(index_data_[beg], index_data_[mid]); 1826 std::swap(value_data_[beg], value_data_[mid]); 1827 } 1828 } else if (len_lef > 0 && len_rig > 0) { 1829 size_type lef_mid, rig_mid; 1830 if (len_lef >= len_rig) { 1831 lef_mid = (beg + mid) / 2; 1832 rig_mid = lower_bound(mid, end, lef_mid); 1833 } else { 1834 rig_mid = (mid + end) / 2; 1835 lef_mid = lower_bound(beg, mid, rig_mid); 1836 } 1837 std::rotate(&index_data_[0] + lef_mid, &index_data_[0] + mid, &index_data_[0] + rig_mid); 1838 std::rotate(&value_data_[0] + lef_mid, &value_data_[0] + mid, &value_data_[0] + rig_mid); 1839 1840 size_type new_mid = lef_mid + rig_mid - mid; 1841 inplace_merge(beg, lef_mid, new_mid); 1842 inplace_merge(new_mid, rig_mid, end); 1843 } 1844 } 1845 1804 1846 // Sorting and summation of duplicates 1805 1847 BOOST_UBLAS_INLINE 1806 1848 void sort () const { … … 1811 1853 const typename array_pair::iterator iunsorted = ipa.begin () + sorted_filled_; 1812 1854 // sort new elements and merge 1813 1855 std::sort (iunsorted, ipa.end ()); 1814 std::inplace_merge (ipa.begin (), iunsorted, ipa.end ()); 1856 // do not use STL algorithm due to compilation problems with recent gcc 1857 inplace_merge(0, sorted_filled_, filled_); 1815 1858 #else 1816 1859 const typename array_pair::iterator iunsorted = ipa.begin (); 1817 1860 std::sort (iunsorted, ipa.end ()); -
boost/numeric/ublas/matrix_sparse.hpp
13 13 #ifndef _BOOST_UBLAS_MATRIX_SPARSE_ 14 14 #define _BOOST_UBLAS_MATRIX_SPARSE_ 15 15 16 #include <typeinfo> 16 17 #include <boost/numeric/ublas/vector_sparse.hpp> 17 18 #include <boost/numeric/ublas/matrix_expression.hpp> 18 19 #include <boost/numeric/ublas/detail/matrix_assign.hpp> … … 54 55 else 55 56 *p = s; 56 57 } 57 58 58 59 public: 59 60 // Construction and destruction 60 61 BOOST_UBLAS_INLINE … … 250 251 * 251 252 * This class represents a matrix by using a \c key to value mapping. The default type is 252 253 * \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 * So, by default a STL map container is used to associate keys and values. The key is computed depending on 254 255 * the layout type \c L as \code key = layout_type::element(i, size1_, j, size2_); \endcode 255 256 * 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 * Limitations: The matrix size must not exceed \f$(size1*size2) < \f$ \code std::limits<std::size_t> \endcode. 257 258 * The \ref find1() and \ref find2() operations have a complexity of at least \f$\mathcal{O}(log(nnz))\f$, depending 258 259 * 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. 260 * Orientation and storage can also be specified, otherwise a row major orientation is used. 261 * It is \b not required by the storage to initialize elements of the matrix. By default, the orientation is \c row_major. 261 262 * 262 263 * \sa fwd.hpp, storage_sparse.hpp 263 264 * … … 427 428 return; 428 429 data ().erase (it); 429 430 } 430 431 431 432 // Zeroing 432 433 BOOST_UBLAS_INLINE 433 434 void clear () { … … 557 558 typedef reverse_iterator_base2<iterator2> reverse_iterator2; 558 559 559 560 // Element lookup 560 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 561 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 561 562 const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const { 562 563 const_subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_))); 563 564 const_subiterator_type it_end (data ().end ()); … … 586 587 } 587 588 return const_iterator1 (*this, rank, i, j, it); 588 589 } 589 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 590 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 590 591 iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) { 591 592 subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_))); 592 593 subiterator_type it_end (data ().end ()); … … 615 616 } 616 617 return iterator1 (*this, rank, i, j, it); 617 618 } 618 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 619 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 619 620 const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const { 620 621 const_subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_))); 621 622 const_subiterator_type it_end (data ().end ()); … … 644 645 } 645 646 return const_iterator2 (*this, rank, i, j, it); 646 647 } 647 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 648 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 648 649 iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) { 649 650 subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_))); 650 651 subiterator_type it_end (data ().end ()); … … 1377 1378 matrix_container<self_type> (), 1378 1379 size1_ (size1), size2_ (size2), data_ () { 1379 1380 data_ [layout_type::size_M (size1_, size2_)] = vector_data_value_type (); 1381 //std::cerr << std::type_info(A).name() << " "; 1382 std::cerr << A::data_value_type(); 1380 1383 } 1381 1384 BOOST_UBLAS_INLINE 1382 1385 mapped_vector_of_mapped_vector (const mapped_vector_of_mapped_vector &m): … … 1509 1512 return; 1510 1513 (*itv).second.erase (it); 1511 1514 } 1512 1515 1513 1516 // Zeroing 1514 1517 BOOST_UBLAS_INLINE 1515 1518 void clear () { … … 1632 1635 subiterator_type it ((*itv).second.find (element2)); 1633 1636 BOOST_UBLAS_CHECK (it != (*itv).second.end (), bad_index ()); 1634 1637 BOOST_UBLAS_CHECK ((*it).first == element2, internal_logic ()); // broken map 1635 1638 1636 1639 return it->second; 1637 1640 } 1638 1641 … … 1647 1650 typedef reverse_iterator_base2<iterator2> reverse_iterator2; 1648 1651 1649 1652 // Element lookup 1650 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 1653 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 1651 1654 const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const { 1652 1655 BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ()); 1653 1656 for (;;) { … … 1662 1665 // advance to the first available major index 1663 1666 size_type M = itv->first; 1664 1667 size_type m; 1665 if (it != it_end) { 1666 m = it->first; 1668 if (it != it_end) { 1669 m = it->first; 1667 1670 } else { 1668 1671 m = layout_type::size_m(size1_, size2_); 1669 1672 } … … 1696 1699 } 1697 1700 } 1698 1701 } 1699 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 1702 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 1700 1703 iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) { 1701 1704 BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ()); 1702 1705 for (;;) { … … 1711 1714 // advance to the first available major index 1712 1715 size_type M = itv->first; 1713 1716 size_type m; 1714 if (it != it_end) { 1715 m = it->first; 1717 if (it != it_end) { 1718 m = it->first; 1716 1719 } else { 1717 1720 m = layout_type::size_m(size1_, size2_); 1718 1721 } … … 1745 1748 } 1746 1749 } 1747 1750 } 1748 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 1751 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 1749 1752 const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const { 1750 1753 BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ()); 1751 1754 for (;;) { … … 1760 1763 // advance to the first available major index 1761 1764 size_type M = itv->first; 1762 1765 size_type m; 1763 if (it != it_end) { 1764 m = it->first; 1766 if (it != it_end) { 1767 m = it->first; 1765 1768 } else { 1766 1769 m = layout_type::size_m(size1_, size2_); 1767 1770 } … … 1794 1797 } 1795 1798 } 1796 1799 } 1797 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 1800 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 1798 1801 iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) { 1799 1802 BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ()); 1800 1803 for (;;) { … … 1809 1812 // advance to the first available major index 1810 1813 size_type M = itv->first; 1811 1814 size_type m; 1812 if (it != it_end) { 1813 m = it->first; 1815 if (it != it_end) { 1816 m = it->first; 1814 1817 } else { 1815 1818 m = layout_type::size_m(size1_, size2_); 1816 1819 } … … 2682 2685 index1_data_ (m.index1_data_), index2_data_ (m.index2_data_), value_data_ (m.value_data_) { 2683 2686 storage_invariants (); 2684 2687 } 2685 2688 2686 2689 BOOST_UBLAS_INLINE 2687 2690 compressed_matrix (const coordinate_matrix<T, L, IB, IA, TA> &m): 2688 2691 matrix_container<self_type> (), … … 2737 2740 size_type nnz () const { 2738 2741 return filled2_; 2739 2742 } 2740 2743 2741 2744 // Storage accessors 2742 2745 BOOST_UBLAS_INLINE 2743 2746 static size_type index_base () { … … 2942 2945 } 2943 2946 storage_invariants (); 2944 2947 } 2945 2948 2946 2949 // Zeroing 2947 2950 BOOST_UBLAS_INLINE 2948 2951 void clear () { … … 3121 3124 typedef reverse_iterator_base2<iterator2> reverse_iterator2; 3122 3125 3123 3126 // Element lookup 3124 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 3127 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 3125 3128 const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const { 3126 3129 for (;;) { 3127 3130 array_size_type address1 (layout_type::index_M (i, j)); … … 3161 3164 } 3162 3165 } 3163 3166 } 3164 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 3167 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 3165 3168 iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) { 3166 3169 for (;;) { 3167 3170 array_size_type address1 (layout_type::index_M (i, j)); … … 3201 3204 } 3202 3205 } 3203 3206 } 3204 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 3207 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 3205 3208 const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const { 3206 3209 for (;;) { 3207 3210 array_size_type address1 (layout_type::index_M (i, j)); … … 3241 3244 } 3242 3245 } 3243 3246 } 3244 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 3247 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 3245 3248 iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) { 3246 3249 for (;;) { 3247 3250 array_size_type address1 (layout_type::index_M (i, j)); … … 3968 3971 BOOST_UBLAS_CHECK (filled2_ <= capacity_, internal_logic ()); 3969 3972 BOOST_UBLAS_CHECK (index1_data_ [filled1_ - 1] == k_based (filled2_), internal_logic ()); 3970 3973 } 3971 3974 3972 3975 size_type size1_; 3973 3976 size_type size2_; 3974 3977 array_size_type capacity_; … … 4208 4211 const_pointer p = find_element (i, j); 4209 4212 if (p) 4210 4213 return *p; 4211 else 4214 else 4212 4215 return zero_; 4213 4216 } 4214 4217 BOOST_UBLAS_INLINE … … 4267 4270 } 4268 4271 storage_invariants (); 4269 4272 } 4270 4273 4271 4274 // Zeroing 4272 4275 BOOST_UBLAS_INLINE 4273 4276 void clear () { … … 4391 4394 m1.swap (m2); 4392 4395 } 4393 4396 4397 // replacement if STL lower bound algorithm for use of inplace_merge 4398 array_size_type lower_bound (array_size_type beg, array_size_type end, array_size_type target) const { 4399 while (end > beg) { 4400 array_size_type mid = (beg + end) / 2; 4401 if (((index1_data_[mid] < index1_data_[target]) || 4402 ((index1_data_[mid] == index1_data_[target]) && 4403 (index2_data_[mid] < index2_data_[target])))) { 4404 beg = mid + 1; 4405 } else { 4406 end = mid; 4407 } 4408 } 4409 return beg; 4410 } 4411 4412 // specialized replacement of STL inplace_merge to avoid compilation 4413 // problems with respect to the array_triple iterator 4414 void inplace_merge (array_size_type beg, array_size_type mid, array_size_type end) const { 4415 array_size_type len_lef = mid - beg; 4416 array_size_type len_rig = end - mid; 4417 4418 if (len_lef == 1 && len_rig == 1) { 4419 if ((index1_data_[mid] < index1_data_[beg]) || 4420 ((index1_data_[mid] == index1_data_[beg]) && (index2_data_[mid] < index2_data_[beg]))) 4421 { 4422 std::swap(index1_data_[beg], index1_data_[mid]); 4423 std::swap(index2_data_[beg], index2_data_[mid]); 4424 std::swap(value_data_[beg], value_data_[mid]); 4425 } 4426 } else if (len_lef > 0 && len_rig > 0) { 4427 array_size_type lef_mid, rig_mid; 4428 if (len_lef >= len_rig) { 4429 lef_mid = (beg + mid) / 2; 4430 rig_mid = lower_bound(mid, end, lef_mid); 4431 } else { 4432 rig_mid = (mid + end) / 2; 4433 lef_mid = lower_bound(beg, mid, rig_mid); 4434 } 4435 std::rotate(&index1_data_[0] + lef_mid, &index1_data_[0] + mid, &index1_data_[0] + rig_mid); 4436 std::rotate(&index2_data_[0] + lef_mid, &index2_data_[0] + mid, &index2_data_[0] + rig_mid); 4437 std::rotate(&value_data_[0] + lef_mid, &value_data_[0] + mid, &value_data_[0] + rig_mid); 4438 4439 array_size_type new_mid = lef_mid + rig_mid - mid; 4440 inplace_merge(beg, lef_mid, new_mid); 4441 inplace_merge(new_mid, rig_mid, end); 4442 } 4443 } 4444 4394 4445 // Sorting and summation of duplicates 4395 4446 BOOST_UBLAS_INLINE 4396 4447 void sort () const { 4397 4448 if (! sorted_ && filled_ > 0) { 4398 4449 typedef index_triple_array<index_array_type, index_array_type, value_array_type> array_triple; 4399 4450 array_triple ita (filled_, index1_data_, index2_data_, value_data_); 4400 #ifndef BOOST_UBLAS_COO_ALWAYS_DO_FULL_SORT4401 4451 const typename array_triple::iterator iunsorted = ita.begin () + sorted_filled_; 4402 4452 // sort new elements and merge 4403 4453 std::sort (iunsorted, ita.end ()); 4404 std::inplace_merge (ita.begin (), iunsorted, ita.end ()); 4405 #else 4406 const typename array_triple::iterator iunsorted = ita.begin (); 4407 std::sort (iunsorted, ita.end ()); 4408 #endif 4454 // do not use STL algorithm due to compilation problems with recent gcc 4455 inplace_merge(0, sorted_filled_, filled_); 4456 4409 4457 // sum duplicates with += and remove 4410 4458 array_size_type filled = 0; 4411 4459 for (array_size_type i = 1; i < filled_; ++ i) { … … 4434 4482 size_type element1 = layout_type::index_M (i, j); 4435 4483 size_type element2 = layout_type::index_m (i, j); 4436 4484 // must maintain sort order 4437 BOOST_UBLAS_CHECK (sorted_ && 4485 BOOST_UBLAS_CHECK (sorted_ && 4438 4486 (filled_ == 0 || 4439 4487 index1_data_ [filled_ - 1] < k_based (element1) || 4440 4488 (index1_data_ [filled_ - 1] == k_based (element1) && index2_data_ [filled_ - 1] < k_based (element2))) … … 4485 4533 typedef reverse_iterator_base2<iterator2> reverse_iterator2; 4486 4534 4487 4535 // Element lookup 4488 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 4536 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 4489 4537 const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const { 4490 4538 sort (); 4491 4539 for (;;) { … … 4526 4574 } 4527 4575 } 4528 4576 } 4529 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 4577 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 4530 4578 iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) { 4531 4579 sort (); 4532 4580 for (;;) { … … 4567 4615 } 4568 4616 } 4569 4617 } 4570 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 4618 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 4571 4619 const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const { 4572 4620 sort (); 4573 4621 for (;;) { … … 4608 4656 } 4609 4657 } 4610 4658 } 4611 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 4659 // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it. 4612 4660 iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) { 4613 4661 sort (); 4614 4662 for (;;) { -
libs/numeric/ublas/test/test_coordinate_vector_inplace_merge.cpp
1 // Copyright (c) 2011 David Bellot 2 // 3 // Distributed under the Boost Software License, Version 1.0. (See 4 // accompanying file LICENSE_1_0.txt or copy at 5 // http://www.boost.org/LICENSE_1_0.txt) 6 7 #ifndef BOOST_UBLAS_NO_ELEMENT_PROXIES 8 # define BOOST_UBLAS_NO_ELEMENT_PROXIES 9 #endif 10 11 #include <boost/numeric/ublas/assignment.hpp> 12 #include <boost/numeric/ublas/vector.hpp> 13 #include <boost/numeric/ublas/vector_sparse.hpp> 14 #include <boost/numeric/ublas/vector_expression.hpp> 15 #include <boost/numeric/ublas/io.hpp> 16 17 #include "libs/numeric/ublas/test/utils.hpp" 18 19 const double TOL = 1e-15; 20 21 22 template <class AE> 23 typename AE::value_type mean_square(const boost::numeric::ublas::vector_expression<AE> &me) { 24 typename AE::value_type s(0); 25 typename AE::size_type i; 26 for (i=0; i!= me().size(); i++) { 27 s += boost::numeric::ublas::scalar_traits<typename AE::value_type>::type_abs(me()(i)); 28 } 29 return s/me().size(); 30 } 31 32 template<typename T> 33 bool check_sortedness(const boost::numeric::ublas::coordinate_vector<T>& vector) { 34 bool result = true; 35 typedef boost::numeric::ublas::coordinate_vector<T> vector_type; 36 typename vector_type::index_array_type idx = vector.index_data(); 37 typename vector_type::size_type size = vector.filled(); 38 39 for (typename vector_type::size_type i = 0; i + 1 < size && result; ++ i) { 40 result &= (idx[i] < idx[i + 1]); 41 } 42 return result; 43 } 44 45 void print_entries(size_t size, 46 const std::vector<size_t>& entries) 47 { 48 std::cerr << "Error entries - Size:" << size << ". Entries: "; 49 for (size_t i = 0; i < entries.size(); ++ i) { 50 std::cerr << entries[i] << "; "; 51 } 52 std::cerr << "\n"; 53 } 54 55 BOOST_UBLAS_TEST_DEF( test_coordinate_vector_inplace_merge_random ) 56 { 57 const size_t max_repeats = 100; 58 const size_t max_size = 100; 59 const size_t dim_var = 10; 60 const size_t nr_entries = 10; 61 62 for (size_t repeats = 1; repeats < max_repeats; ++repeats ) { 63 for (size_t size = 1; size < max_size; size += 5) { 64 size_t size_vec = size + rand() % dim_var; 65 66 boost::numeric::ublas::coordinate_vector<double> vector_coord(size_vec); 67 boost::numeric::ublas::vector<double> vector_dense(size_vec, 0); 68 69 vector_coord.sort(); 70 71 std::vector<size_t> entries; 72 for (int entry = 0; entry < nr_entries; ++ entry) { 73 int x = rand() % size_vec; 74 entries.push_back(x); 75 vector_coord.append_element(x, 1); 76 vector_dense(x) += 1; 77 } 78 vector_coord.sort(); 79 80 { 81 bool sorted = check_sortedness(vector_coord); 82 bool identical = mean_square(vector_coord - vector_dense) < TOL; 83 if (!(sorted && identical)) { 84 print_entries(size_vec, entries); 85 } 86 BOOST_UBLAS_TEST_CHECK( check_sortedness(vector_coord) ); 87 BOOST_UBLAS_TEST_CHECK( mean_square(vector_coord - vector_dense) < TOL); 88 } 89 90 for (int entry = 0; entry < nr_entries; ++ entry) { 91 int x = rand() % size_vec; 92 entries.push_back(x); 93 vector_coord(x) += 1; 94 vector_dense(x) += 1; 95 vector_coord.sort(); 96 } 97 98 { 99 bool sorted = check_sortedness(vector_coord); 100 bool identical = mean_square(vector_coord - vector_dense) < TOL; 101 if (!(sorted && identical)) { 102 print_entries(size_vec, entries); 103 } 104 BOOST_UBLAS_TEST_CHECK( sorted ); 105 BOOST_UBLAS_TEST_CHECK( identical ); 106 } 107 } 108 } 109 } 110 111 int main() 112 { 113 BOOST_UBLAS_TEST_BEGIN(); 114 115 BOOST_UBLAS_TEST_DO( test_coordinate_vector_inplace_merge_random ); 116 117 BOOST_UBLAS_TEST_END(); 118 119 return EXIT_SUCCESS;; 120 } -
libs/numeric/ublas/test/test_coordinate_matrix_inplace_merge.cpp
1 // Copyright (c) 2011 David Bellot 2 // 3 // Distributed under the Boost Software License, Version 1.0. (See 4 // accompanying file LICENSE_1_0.txt or copy at 5 // http://www.boost.org/LICENSE_1_0.txt) 6 7 #ifndef BOOST_UBLAS_NO_ELEMENT_PROXIES 8 # define BOOST_UBLAS_NO_ELEMENT_PROXIES 9 #endif 10 11 #include <boost/numeric/ublas/assignment.hpp> 12 #include <boost/numeric/ublas/matrix.hpp> 13 #include <boost/numeric/ublas/matrix_sparse.hpp> 14 #include <boost/numeric/ublas/matrix_expression.hpp> 15 #include <boost/numeric/ublas/io.hpp> 16 17 #include "libs/numeric/ublas/test/utils.hpp" 18 19 using std::cout; 20 using std::endl; 21 22 const double TOL = 1e-15; 23 24 template <class AE> 25 typename AE::value_type mean_square(const boost::numeric::ublas::matrix_expression<AE> &me) { 26 typename AE::value_type s(0); 27 typename AE::size_type i, j; 28 for (i=0; i!= me().size1(); i++) { 29 for (j=0; j!= me().size2(); j++) { 30 s += boost::numeric::ublas::scalar_traits<typename AE::value_type>::type_abs(me()(i,j)); 31 } 32 } 33 return s/me().size1()*me().size2(); 34 } 35 36 template<typename T> 37 bool check_sortedness(const boost::numeric::ublas::coordinate_matrix<T>& matrix) { 38 bool result = true; 39 typedef boost::numeric::ublas::coordinate_matrix<T> matrix_type; 40 typename matrix_type::index_array_type i1 = matrix.index1_data(); 41 typename matrix_type::index_array_type i2 = matrix.index2_data(); 42 typename matrix_type::array_size_type size = matrix.filled(); 43 44 for (typename matrix_type::array_size_type i = 0; i + 1 < size && result; ++ i) { 45 result &= ( (i1[i] < i1[i + 1]) || 46 ((i1[i] == i1[i]) && 47 (i2[i] < i2[i + 1])) ); 48 49 } 50 return result; 51 } 52 53 void print_entries(size_t size_x, size_t size_y, 54 const std::vector<std::pair<size_t, size_t> >& entries) 55 { 56 std::cerr << "Error - Size:" << size_x << " x " << size_y << ". Entries: "; 57 for (size_t i = 0; i < entries.size(); ++ i) { 58 std::cerr << entries[i].first << ", " << entries[i].second << "; "; 59 } 60 std::cerr << "\n"; 61 } 62 63 64 BOOST_UBLAS_TEST_DEF( test_coordinate_matrix_inplace_merge_random ) 65 { 66 const size_t max_repeats = 100; 67 const size_t max_size = 100; 68 const size_t dim_var = 10; 69 const size_t nr_entries = 10; 70 71 for (size_t repeats = 1; repeats < max_repeats; ++repeats ) { 72 for (size_t size = 1; size < max_size; size += 5) { 73 size_t size_x = size + rand() % dim_var; 74 size_t size_y = size + rand() % dim_var; 75 76 boost::numeric::ublas::coordinate_matrix<double> matrix_coord(size_x, size_y); 77 boost::numeric::ublas::matrix<double> matrix_dense(size_x, size_y, 0); 78 79 matrix_coord.sort(); 80 81 std::vector<std::pair<size_t, size_t> > entries; 82 for (int entry = 0; entry < nr_entries; ++ entry) { 83 int x = rand() % size_x; 84 int y = rand() % size_y; 85 entries.push_back(std::make_pair(x, y)); 86 matrix_coord.append_element(x, y, 1); 87 matrix_dense(x, y) += 1; 88 } 89 matrix_coord.sort(); 90 91 { 92 bool sorted = check_sortedness(matrix_coord); 93 bool identical = mean_square(matrix_coord - matrix_dense) < TOL; 94 if (!(sorted && identical)) { 95 print_entries(size_x, size_y, entries); 96 } 97 BOOST_UBLAS_TEST_CHECK( check_sortedness(matrix_coord) ); 98 BOOST_UBLAS_TEST_CHECK( mean_square(matrix_coord - matrix_dense) < TOL); 99 } 100 101 for (int entry = 0; entry < nr_entries; ++ entry) { 102 int x = rand() % size_x; 103 int y = rand() % size_y; 104 entries.push_back(std::make_pair(x, y)); 105 matrix_coord(x, y) += 1; 106 matrix_dense(x, y) += 1; 107 matrix_coord.sort(); 108 } 109 110 { 111 bool sorted = check_sortedness(matrix_coord); 112 bool identical = mean_square(matrix_coord - matrix_dense) < TOL; 113 if (!(sorted && identical)) { 114 print_entries(size_x, size_y, entries); 115 } 116 BOOST_UBLAS_TEST_CHECK( sorted ); 117 BOOST_UBLAS_TEST_CHECK( identical ); 118 } 119 } 120 } 121 } 122 123 int main() 124 { 125 BOOST_UBLAS_TEST_BEGIN(); 126 127 BOOST_UBLAS_TEST_DO( test_coordinate_matrix_inplace_merge_random ); 128 129 BOOST_UBLAS_TEST_END(); 130 131 return EXIT_SUCCESS;; 132 }
