Ticket #7363: custom_inplace_merge.diff

File custom_inplace_merge.diff, 33.5 KB (added by anonymous, 10 years ago)

  • boost/numeric/ublas/vector_sparse.hpp

     
    262262
    263263    /** \brief Index map based sparse vector
    264264     *
    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
    266266     * \c std::map<size_t, T> or \c map_array<size_t, T>. This means that only non-zero elements
    267267     * are effectively stored.
    268268     *
    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
    271271     * \f$k = v_{i_1}\f$ and \f$k + 1 = v_{i_2}\f$ of the container, holds \f$i_1 < i_2\f$.
    272272     *
    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
    274274     * \c map_std<std::size_t, T>. The latter is equivalent to \c std::map<std::size_t, T>.
    275275     *
    276276     * \tparam T the type of object stored in the vector (like double, float, complex, etc...)
    … …  
    771771
    772772
    773773    // Thanks to Kresimir Fresl for extending this to cover different index bases.
    774    
     774
    775775    /** \brief Compressed array based sparse vector
    776776     *
    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
    779779     * and there is at most one entry for each index. Inserting an element can be time consuming.
    780780     * If the vector contains a few zero entries, then it is better to have a normal vector.
    781781     * If the vector has a very high dimension with a few non-zero values, then this vector is
    782782     * very memory efficient (at the cost of a few more computations).
    783783     *
    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
    786786     * 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$.
    787787     *
    788788     * Supported parameters for the adapted array (indices and values) are \c unbounded_array<> ,
    … …  
    14141414
    14151415    /** \brief Coordimate array based sparse vector
    14161416     *
    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
    14201420     * 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
    14221422     * be longer in specific cases than a \c compressed_vector.
    14231423     *
    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
    14261426     * 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$.
    14271427     *
    14281428     * Supported parameters for the adapted array (indices and values) are \c unbounded_array<> ,
    … …  
    18011801            v1.swap (v2);
    18021802        }
    18031803
     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
    18041846        // Sorting and summation of duplicates
    18051847        BOOST_UBLAS_INLINE
    18061848        void sort () const {
    … …  
    18101852                const typename array_pair::iterator iunsorted = ipa.begin () + sorted_filled_;
    18111853                // sort new elements and merge
    18121854                std::sort (iunsorted, ipa.end ());
    1813                 std::inplace_merge (ipa.begin (), iunsorted, ipa.end ());
     1855                // do not use STL algorithm due to compilation problems with recent gcc
     1856                inplace_merge(0, sorted_filled_, filled_);
    18141857
    18151858                // sum duplicates with += and remove
    18161859                size_type filled = 0;

  • boost/numeric/ublas/matrix_sparse.hpp

     
    1313#ifndef _BOOST_UBLAS_MATRIX_SPARSE_
    1414#define _BOOST_UBLAS_MATRIX_SPARSE_
    1515
     16#include <typeinfo>
    1617#include <boost/numeric/ublas/vector_sparse.hpp>
    1718#include <boost/numeric/ublas/matrix_expression.hpp>
    1819#include <boost/numeric/ublas/detail/matrix_assign.hpp>
     20#include <boost/numeric/ublas/inplace_merge.hpp>
    1921#if BOOST_UBLAS_TYPE_CHECK
    2022#include <boost/numeric/ublas/matrix.hpp>
    2123#endif
    … …  
    5456            else
    5557                *p = s;
    5658        }
    57        
     59
    5860    public:
    5961        // Construction and destruction
    6062        BOOST_UBLAS_INLINE
    … …  
    250252    *
    251253    * This class represents a matrix by using a \c key to value mapping. The default type is
    252254    * \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 
     255    * So, by default a STL map container is used to associate keys and values. The key is computed depending on
    254256    * the layout type \c L as \code key = layout_type::element(i, size1_, j, size2_); \endcode
    255257    * 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. 
     258    * Limitations: The matrix size must not exceed \f$(size1*size2) < \f$ \code std::limits<std::size_t> \endcode.
    257259    * The \ref find1() and \ref find2() operations have a complexity of at least \f$\mathcal{O}(log(nnz))\f$, depending
    258260    * on the efficiency of \c std::lower_bound on the key set of the map.
    259     * Orientation and storage can also be specified, otherwise a row major orientation is used. 
    260     * It is \b not required by the storage to initialize elements of the matrix. By default, the orientation is \c row_major. 
     261    * Orientation and storage can also be specified, otherwise a row major orientation is used.
     262    * It is \b not required by the storage to initialize elements of the matrix. By default, the orientation is \c row_major.
    261263    *
    262264    * \sa fwd.hpp, storage_sparse.hpp
    263265    *
    … …  
    427429                return;
    428430            data ().erase (it);
    429431        }
    430        
     432
    431433        // Zeroing
    432434        BOOST_UBLAS_INLINE
    433435        void clear () {
    … …  
    557559        typedef reverse_iterator_base2<iterator2> reverse_iterator2;
    558560
    559561        // Element lookup
    560         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     562        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    561563        const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const {
    562564            const_subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_)));
    563565            const_subiterator_type it_end (data ().end ());
    … …  
    586588            }
    587589            return const_iterator1 (*this, rank, i, j, it);
    588590        }
    589         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     591        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    590592        iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) {
    591593            subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_)));
    592594            subiterator_type it_end (data ().end ());
    … …  
    615617            }
    616618            return iterator1 (*this, rank, i, j, it);
    617619        }
    618         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     620        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    619621        const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const {
    620622            const_subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_)));
    621623            const_subiterator_type it_end (data ().end ());
    … …  
    644646            }
    645647            return const_iterator2 (*this, rank, i, j, it);
    646648        }
    647         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     649        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    648650        iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) {
    649651            subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_)));
    650652            subiterator_type it_end (data ().end ());
    … …  
    13771379            matrix_container<self_type> (),
    13781380            size1_ (size1), size2_ (size2), data_ () {
    13791381            data_ [layout_type::size_M (size1_, size2_)] = vector_data_value_type ();
     1382            //std::cerr << std::type_info(A).name() << " ";
     1383            std::cerr << A::data_value_type();
    13801384        }
    13811385        BOOST_UBLAS_INLINE
    13821386        mapped_vector_of_mapped_vector (const mapped_vector_of_mapped_vector &m):
    … …  
    15091513                return;
    15101514            (*itv).second.erase (it);
    15111515        }
    1512        
     1516
    15131517        // Zeroing
    15141518        BOOST_UBLAS_INLINE
    15151519        void clear () {
    … …  
    16321636            subiterator_type it ((*itv).second.find (element2));
    16331637            BOOST_UBLAS_CHECK (it != (*itv).second.end (), bad_index ());
    16341638            BOOST_UBLAS_CHECK ((*it).first == element2, internal_logic ());    // broken map
    1635            
     1639
    16361640            return it->second;
    16371641        }
    16381642
    … …  
    16471651        typedef reverse_iterator_base2<iterator2> reverse_iterator2;
    16481652
    16491653        // Element lookup
    1650         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     1654        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    16511655        const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const {
    16521656            BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ());
    16531657            for (;;) {
    … …  
    16621666                    // advance to the first available major index
    16631667                    size_type M = itv->first;
    16641668                    size_type m;
    1665                     if (it != it_end) { 
    1666                         m = it->first; 
     1669                    if (it != it_end) {
     1670                        m = it->first;
    16671671                    } else {
    16681672                        m = layout_type::size_m(size1_, size2_);
    16691673                    }
    … …  
    16961700                }
    16971701            }
    16981702        }
    1699         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     1703        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    17001704        iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) {
    17011705            BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ());
    17021706            for (;;) {
    … …  
    17111715                    // advance to the first available major index
    17121716                    size_type M = itv->first;
    17131717                    size_type m;
    1714                     if (it != it_end) { 
    1715                         m = it->first; 
     1718                    if (it != it_end) {
     1719                        m = it->first;
    17161720                    } else {
    17171721                        m = layout_type::size_m(size1_, size2_);
    17181722                    }
    … …  
    17451749                }
    17461750            }
    17471751        }
    1748         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     1752        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    17491753        const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const {
    17501754            BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ());
    17511755            for (;;) {
    … …  
    17601764                    // advance to the first available major index
    17611765                    size_type M = itv->first;
    17621766                    size_type m;
    1763                     if (it != it_end) { 
    1764                         m = it->first; 
     1767                    if (it != it_end) {
     1768                        m = it->first;
    17651769                    } else {
    17661770                        m = layout_type::size_m(size1_, size2_);
    17671771                    }
    … …  
    17941798                }
    17951799            }
    17961800        }
    1797         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     1801        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    17981802        iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) {
    17991803            BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ());
    18001804            for (;;) {
    … …  
    18091813                    // advance to the first available major index
    18101814                    size_type M = itv->first;
    18111815                    size_type m;
    1812                     if (it != it_end) { 
    1813                         m = it->first; 
     1816                    if (it != it_end) {
     1817                        m = it->first;
    18141818                    } else {
    18151819                        m = layout_type::size_m(size1_, size2_);
    18161820                    }
    … …  
    26822686            index1_data_ (m.index1_data_), index2_data_ (m.index2_data_), value_data_ (m.value_data_) {
    26832687            storage_invariants ();
    26842688        }
    2685          
     2689
    26862690        BOOST_UBLAS_INLINE
    26872691        compressed_matrix (const coordinate_matrix<T, L, IB, IA, TA> &m):
    26882692            matrix_container<self_type> (),
    … …  
    27372741        size_type nnz () const {
    27382742            return filled2_;
    27392743        }
    2740        
     2744
    27412745        // Storage accessors
    27422746        BOOST_UBLAS_INLINE
    27432747        static size_type index_base () {
    … …  
    29422946            }
    29432947            storage_invariants ();
    29442948        }
    2945        
     2949
    29462950        // Zeroing
    29472951        BOOST_UBLAS_INLINE
    29482952        void clear () {
    … …  
    31213125        typedef reverse_iterator_base2<iterator2> reverse_iterator2;
    31223126
    31233127        // Element lookup
    3124         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     3128        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    31253129        const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const {
    31263130            for (;;) {
    31273131                array_size_type address1 (layout_type::index_M (i, j));
    … …  
    31613165                }
    31623166            }
    31633167        }
    3164         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     3168        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    31653169        iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) {
    31663170            for (;;) {
    31673171                array_size_type address1 (layout_type::index_M (i, j));
    … …  
    32013205                }
    32023206            }
    32033207        }
    3204         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     3208        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    32053209        const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const {
    32063210            for (;;) {
    32073211                array_size_type address1 (layout_type::index_M (i, j));
    … …  
    32413245                }
    32423246            }
    32433247        }
    3244         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     3248        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    32453249        iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) {
    32463250            for (;;) {
    32473251                array_size_type address1 (layout_type::index_M (i, j));
    … …  
    39683972            BOOST_UBLAS_CHECK (filled2_ <= capacity_, internal_logic ());
    39693973            BOOST_UBLAS_CHECK (index1_data_ [filled1_ - 1] == k_based (filled2_), internal_logic ());
    39703974        }
    3971        
     3975
    39723976        size_type size1_;
    39733977        size_type size2_;
    39743978        array_size_type capacity_;
    … …  
    42084212            const_pointer p = find_element (i, j);
    42094213            if (p)
    42104214                return *p;
    4211             else 
     4215            else
    42124216                return zero_;
    42134217        }
    42144218        BOOST_UBLAS_INLINE
    … …  
    42674271            }
    42684272            storage_invariants ();
    42694273        }
    4270        
     4274
    42714275        // Zeroing
    42724276        BOOST_UBLAS_INLINE
    42734277        void clear () {
    … …  
    43914395            m1.swap (m2);
    43924396        }
    43934397
     4398        // replacement if STL lower bound algorithm for use of inplace_merge
     4399        array_size_type lower_bound (array_size_type beg, array_size_type end, array_size_type target) const {
     4400          while (end > beg) {
     4401            array_size_type mid = (beg + end) / 2;
     4402            if (((index1_data_[mid] < index1_data_[target]) ||
     4403                 ((index1_data_[mid] == index1_data_[target]) &&
     4404                  (index2_data_[mid] < index2_data_[target])))) {
     4405              beg = mid + 1;
     4406            } else {
     4407              end = mid;
     4408            }
     4409          }
     4410          return beg;
     4411        }
     4412
     4413        // specialized replacement of STL inplace_merge to avoid compilation
     4414        // problems with respect to the array_triple iterator
     4415        void inplace_merge (array_size_type beg, array_size_type mid, array_size_type end) const {
     4416          array_size_type len_lef = mid - beg;
     4417          array_size_type len_rig = end - mid;
     4418
     4419          if (len_lef == 1 && len_rig == 1) {
     4420            if ((index1_data_[mid] < index1_data_[beg]) ||
     4421                ((index1_data_[mid] == index1_data_[beg]) && (index2_data_[mid] < index2_data_[beg])))
     4422              {
     4423                std::swap(index1_data_[beg], index1_data_[mid]);
     4424                std::swap(index2_data_[beg], index2_data_[mid]);
     4425                std::swap(value_data_[beg], value_data_[mid]);
     4426              }
     4427          } else if (len_lef > 0 && len_rig > 0) {
     4428            array_size_type lef_mid, rig_mid;
     4429            if (len_lef >= len_rig) {
     4430              lef_mid = (beg + mid) / 2;
     4431              rig_mid = lower_bound(mid, end, lef_mid);
     4432            } else {
     4433              rig_mid = (mid + end) / 2;
     4434              lef_mid = lower_bound(beg, mid, rig_mid);
     4435            }
     4436            std::rotate(&index1_data_[0] + lef_mid, &index1_data_[0] + mid, &index1_data_[0] + rig_mid);
     4437            std::rotate(&index2_data_[0] + lef_mid, &index2_data_[0] + mid, &index2_data_[0] + rig_mid);
     4438            std::rotate(&value_data_[0] + lef_mid, &value_data_[0] + mid, &value_data_[0] + rig_mid);
     4439
     4440            array_size_type new_mid = lef_mid + rig_mid - mid;
     4441            inplace_merge(beg, lef_mid, new_mid);
     4442            inplace_merge(new_mid, rig_mid, end);
     4443          }
     4444        }
     4445
    43944446        // Sorting and summation of duplicates
    43954447        BOOST_UBLAS_INLINE
    43964448        void sort () const {
    43974449            if (! sorted_ && filled_ > 0) {
    43984450                typedef index_triple_array<index_array_type, index_array_type, value_array_type> array_triple;
    43994451                array_triple ita (filled_, index1_data_, index2_data_, value_data_);
    4400 #ifndef BOOST_UBLAS_COO_ALWAYS_DO_FULL_SORT
    44014452                const typename array_triple::iterator iunsorted = ita.begin () + sorted_filled_;
    44024453                // sort new elements and merge
    44034454                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               
     4455                // do not use STL algorithm due to compilation problems with recent gcc
     4456                inplace_merge(0, sorted_filled_, filled_);
     4457
    44094458                // sum duplicates with += and remove
    44104459                array_size_type filled = 0;
    44114460                for (array_size_type i = 1; i < filled_; ++ i) {
    … …  
    44344483            size_type element1 = layout_type::index_M (i, j);
    44354484            size_type element2 = layout_type::index_m (i, j);
    44364485            // must maintain sort order
    4437             BOOST_UBLAS_CHECK (sorted_ && 
     4486            BOOST_UBLAS_CHECK (sorted_ &&
    44384487                    (filled_ == 0 ||
    44394488                    index1_data_ [filled_ - 1] < k_based (element1) ||
    44404489                    (index1_data_ [filled_ - 1] == k_based (element1) && index2_data_ [filled_ - 1] < k_based (element2)))
    … …  
    44854534        typedef reverse_iterator_base2<iterator2> reverse_iterator2;
    44864535
    44874536        // Element lookup
    4488         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     4537        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    44894538        const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const {
    44904539            sort ();
    44914540            for (;;) {
    … …  
    45264575                }
    45274576            }
    45284577        }
    4529         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     4578        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    45304579        iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) {
    45314580            sort ();
    45324581            for (;;) {
    … …  
    45674616                }
    45684617            }
    45694618        }
    4570         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     4619        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    45714620        const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const {
    45724621            sort ();
    45734622            for (;;) {
    … …  
    46084657                }
    46094658            }
    46104659        }
    4611         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     4660        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    46124661        iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) {
    46134662            sort ();
    46144663            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
     19const double TOL = 1e-15;
     20
     21
     22template <class AE>
     23typename 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
     32template<typename T>
     33bool 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
     45void 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
     55BOOST_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
     111int 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
     19using std::cout;
     20using std::endl;
     21
     22const double TOL = 1e-15;
     23
     24template <class AE>
     25typename 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
     36template<typename T>
     37bool 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
     53void 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
     64BOOST_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
     123int 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}