Ticket #7363: custom_inplace_merge_r80623.diff

File custom_inplace_merge_r80623.diff, 32.4 KB (added by j.ungermann@…, 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 {
    … …  
    18111853                const typename array_pair::iterator iunsorted = ipa.begin () + sorted_filled_;
    18121854                // sort new elements and merge
    18131855                std::sort (iunsorted, ipa.end ());
    1814                 std::inplace_merge (ipa.begin (), iunsorted, ipa.end ());
     1856                inplace_merge(0, sorted_filled_, filled_);
    18151857#else
    18161858                const typename array_pair::iterator iunsorted = ipa.begin ();
    18171859                std::sort (iunsorted, ipa.end ());
    1818 #endif               
     1860#endif
    18191861
    18201862                // sum duplicates with += and remove
    18211863                size_type filled = 0;

  • boost/numeric/ublas/matrix_sparse.hpp

     
    5454            else
    5555                *p = s;
    5656        }
    57        
     57
    5858    public:
    5959        // Construction and destruction
    6060        BOOST_UBLAS_INLINE
    … …  
    250250    *
    251251    * This class represents a matrix by using a \c key to value mapping. The default type is
    252252    * \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 
     253    * So, by default a STL map container is used to associate keys and values. The key is computed depending on
    254254    * the layout type \c L as \code key = layout_type::element(i, size1_, j, size2_); \endcode
    255255    * 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. 
     256    * Limitations: The matrix size must not exceed \f$(size1*size2) < \f$ \code std::limits<std::size_t> \endcode.
    257257    * The \ref find1() and \ref find2() operations have a complexity of at least \f$\mathcal{O}(log(nnz))\f$, depending
    258258    * 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. 
     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.
    261261    *
    262262    * \sa fwd.hpp, storage_sparse.hpp
    263263    *
    … …  
    427427                return;
    428428            data ().erase (it);
    429429        }
    430        
     430
    431431        // Zeroing
    432432        BOOST_UBLAS_INLINE
    433433        void clear () {
    … …  
    557557        typedef reverse_iterator_base2<iterator2> reverse_iterator2;
    558558
    559559        // Element lookup
    560         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     560        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    561561        const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const {
    562562            const_subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_)));
    563563            const_subiterator_type it_end (data ().end ());
    … …  
    586586            }
    587587            return const_iterator1 (*this, rank, i, j, it);
    588588        }
    589         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     589        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    590590        iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) {
    591591            subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_)));
    592592            subiterator_type it_end (data ().end ());
    … …  
    615615            }
    616616            return iterator1 (*this, rank, i, j, it);
    617617        }
    618         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     618        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    619619        const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const {
    620620            const_subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_)));
    621621            const_subiterator_type it_end (data ().end ());
    … …  
    644644            }
    645645            return const_iterator2 (*this, rank, i, j, it);
    646646        }
    647         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     647        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    648648        iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) {
    649649            subiterator_type it (data ().lower_bound (layout_type::address (i, size1_, j, size2_)));
    650650            subiterator_type it_end (data ().end ());
    … …  
    15091509                return;
    15101510            (*itv).second.erase (it);
    15111511        }
    1512        
     1512
    15131513        // Zeroing
    15141514        BOOST_UBLAS_INLINE
    15151515        void clear () {
    … …  
    16321632            subiterator_type it ((*itv).second.find (element2));
    16331633            BOOST_UBLAS_CHECK (it != (*itv).second.end (), bad_index ());
    16341634            BOOST_UBLAS_CHECK ((*it).first == element2, internal_logic ());    // broken map
    1635            
     1635
    16361636            return it->second;
    16371637        }
    16381638
    … …  
    16471647        typedef reverse_iterator_base2<iterator2> reverse_iterator2;
    16481648
    16491649        // Element lookup
    1650         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     1650        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    16511651        const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const {
    16521652            BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ());
    16531653            for (;;) {
    … …  
    16621662                    // advance to the first available major index
    16631663                    size_type M = itv->first;
    16641664                    size_type m;
    1665                     if (it != it_end) { 
    1666                         m = it->first; 
     1665                    if (it != it_end) {
     1666                        m = it->first;
    16671667                    } else {
    16681668                        m = layout_type::size_m(size1_, size2_);
    16691669                    }
    … …  
    16961696                }
    16971697            }
    16981698        }
    1699         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     1699        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    17001700        iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) {
    17011701            BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ());
    17021702            for (;;) {
    … …  
    17111711                    // advance to the first available major index
    17121712                    size_type M = itv->first;
    17131713                    size_type m;
    1714                     if (it != it_end) { 
    1715                         m = it->first; 
     1714                    if (it != it_end) {
     1715                        m = it->first;
    17161716                    } else {
    17171717                        m = layout_type::size_m(size1_, size2_);
    17181718                    }
    … …  
    17451745                }
    17461746            }
    17471747        }
    1748         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     1748        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    17491749        const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const {
    17501750            BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ());
    17511751            for (;;) {
    … …  
    17601760                    // advance to the first available major index
    17611761                    size_type M = itv->first;
    17621762                    size_type m;
    1763                     if (it != it_end) { 
    1764                         m = it->first; 
     1763                    if (it != it_end) {
     1764                        m = it->first;
    17651765                    } else {
    17661766                        m = layout_type::size_m(size1_, size2_);
    17671767                    }
    … …  
    17941794                }
    17951795            }
    17961796        }
    1797         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     1797        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    17981798        iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) {
    17991799            BOOST_UBLAS_CHECK (data ().begin () != data ().end (), internal_logic ());
    18001800            for (;;) {
    … …  
    18091809                    // advance to the first available major index
    18101810                    size_type M = itv->first;
    18111811                    size_type m;
    1812                     if (it != it_end) { 
    1813                         m = it->first; 
     1812                    if (it != it_end) {
     1813                        m = it->first;
    18141814                    } else {
    18151815                        m = layout_type::size_m(size1_, size2_);
    18161816                    }
    … …  
    26822682            index1_data_ (m.index1_data_), index2_data_ (m.index2_data_), value_data_ (m.value_data_) {
    26832683            storage_invariants ();
    26842684        }
    2685          
     2685
    26862686        BOOST_UBLAS_INLINE
    26872687        compressed_matrix (const coordinate_matrix<T, L, IB, IA, TA> &m):
    26882688            matrix_container<self_type> (),
    … …  
    27372737        size_type nnz () const {
    27382738            return filled2_;
    27392739        }
    2740        
     2740
    27412741        // Storage accessors
    27422742        BOOST_UBLAS_INLINE
    27432743        static size_type index_base () {
    … …  
    29422942            }
    29432943            storage_invariants ();
    29442944        }
    2945        
     2945
    29462946        // Zeroing
    29472947        BOOST_UBLAS_INLINE
    29482948        void clear () {
    … …  
    31213121        typedef reverse_iterator_base2<iterator2> reverse_iterator2;
    31223122
    31233123        // Element lookup
    3124         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     3124        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    31253125        const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const {
    31263126            for (;;) {
    31273127                array_size_type address1 (layout_type::index_M (i, j));
    … …  
    31613161                }
    31623162            }
    31633163        }
    3164         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     3164        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    31653165        iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) {
    31663166            for (;;) {
    31673167                array_size_type address1 (layout_type::index_M (i, j));
    … …  
    32013201                }
    32023202            }
    32033203        }
    3204         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     3204        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    32053205        const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const {
    32063206            for (;;) {
    32073207                array_size_type address1 (layout_type::index_M (i, j));
    … …  
    32413241                }
    32423242            }
    32433243        }
    3244         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
     3244        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.
    32453245        iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) {
    32463246            for (;;) {
    32473247                array_size_type address1 (layout_type::index_M (i, j));
    … …  
    39683968            BOOST_UBLAS_CHECK (filled2_ <= capacity_, internal_logic ());
    39693969            BOOST_UBLAS_CHECK (index1_data_ [filled1_ - 1] == k_based (filled2_), internal_logic ());
    39703970        }
    3971        
     3971
    39723972        size_type size1_;
    39733973        size_type size2_;
    39743974        array_size_type capacity_;
    … …  
    42084208            const_pointer p = find_element (i, j);
    42094209            if (p)
    42104210                return *p;
    4211             else 
     4211            else
    42124212                return zero_;
    42134213        }
    42144214        BOOST_UBLAS_INLINE
    … …  
    42674267            }
    42684268            storage_invariants ();
    42694269        }
    4270        
     4270
    42714271        // Zeroing
    42724272        BOOST_UBLAS_INLINE
    42734273        void clear () {
    … …  
    43914391            m1.swap (m2);
    43924392        }
    43934393
     4394        // replacement if STL lower bound algorithm for use of inplace_merge
     4395        array_size_type lower_bound (array_size_type beg, array_size_type end, array_size_type target) const {
     4396          while (end > beg) {
     4397            array_size_type mid = (beg + end) / 2;
     4398            if (((index1_data_[mid] < index1_data_[target]) ||
     4399                 ((index1_data_[mid] == index1_data_[target]) &&
     4400                  (index2_data_[mid] < index2_data_[target])))) {
     4401              beg = mid + 1;
     4402            } else {
     4403              end = mid;
     4404            }
     4405          }
     4406          return beg;
     4407        }
     4408
     4409        // specialized replacement of STL inplace_merge to avoid compilation
     4410        // problems with respect to the array_triple iterator
     4411        void inplace_merge (array_size_type beg, array_size_type mid, array_size_type end) const {
     4412          array_size_type len_lef = mid - beg;
     4413          array_size_type len_rig = end - mid;
     4414
     4415          if (len_lef == 1 && len_rig == 1) {
     4416            if ((index1_data_[mid] < index1_data_[beg]) ||
     4417                ((index1_data_[mid] == index1_data_[beg]) && (index2_data_[mid] < index2_data_[beg])))
     4418              {
     4419                std::swap(index1_data_[beg], index1_data_[mid]);
     4420                std::swap(index2_data_[beg], index2_data_[mid]);
     4421                std::swap(value_data_[beg], value_data_[mid]);
     4422              }
     4423          } else if (len_lef > 0 && len_rig > 0) {
     4424            array_size_type lef_mid, rig_mid;
     4425            if (len_lef >= len_rig) {
     4426              lef_mid = (beg + mid) / 2;
     4427              rig_mid = lower_bound(mid, end, lef_mid);
     4428            } else {
     4429              rig_mid = (mid + end) / 2;
     4430              lef_mid = lower_bound(beg, mid, rig_mid);
     4431            }
     4432            std::rotate(&index1_data_[0] + lef_mid, &index1_data_[0] + mid, &index1_data_[0] + rig_mid);
     4433            std::rotate(&index2_data_[0] + lef_mid, &index2_data_[0] + mid, &index2_data_[0] + rig_mid);
     4434            std::rotate(&value_data_[0] + lef_mid, &value_data_[0] + mid, &value_data_[0] + rig_mid);
     4435
     4436            array_size_type new_mid = lef_mid + rig_mid - mid;
     4437            inplace_merge(beg, lef_mid, new_mid);
     4438            inplace_merge(new_mid, rig_mid, end);
     4439          }
     4440        }
     4441
    43944442        // Sorting and summation of duplicates
    43954443        BOOST_UBLAS_INLINE
    43964444        void sort () const {
    … …  
    44014449                const typename array_triple::iterator iunsorted = ita.begin () + sorted_filled_;
    44024450                // sort new elements and merge
    44034451                std::sort (iunsorted, ita.end ());
    4404                 std::inplace_merge (ita.begin (), iunsorted, ita.end ());
     4452                inplace_merge(0, sorted_filled_, filled_);
    44054453#else
    44064454                const typename array_triple::iterator iunsorted = ita.begin ();
    44074455                std::sort (iunsorted, ita.end ());
    4408 #endif               
     4456#endif
    44094457                // sum duplicates with += and remove
    44104458                array_size_type filled = 0;
    44114459                for (array_size_type i = 1; i < filled_; ++ i) {
    … …  
    44344482            size_type element1 = layout_type::index_M (i, j);
    44354483            size_type element2 = layout_type::index_m (i, j);
    44364484            // must maintain sort order
    4437             BOOST_UBLAS_CHECK (sorted_ && 
     4485            BOOST_UBLAS_CHECK (sorted_ &&
    44384486                    (filled_ == 0 ||
    44394487                    index1_data_ [filled_ - 1] < k_based (element1) ||
    44404488                    (index1_data_ [filled_ - 1] == k_based (element1) && index2_data_ [filled_ - 1] < k_based (element2)))
    … …  
    44854533        typedef reverse_iterator_base2<iterator2> reverse_iterator2;
    44864534
    44874535        // 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.
    44894537        const_iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) const {
    44904538            sort ();
    44914539            for (;;) {
    … …  
    45264574                }
    45274575            }
    45284576        }
    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.
    45304578        iterator1 find1 (int rank, size_type i, size_type j, int direction = 1) {
    45314579            sort ();
    45324580            for (;;) {
    … …  
    45674615                }
    45684616            }
    45694617        }
    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.
    45714619        const_iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) const {
    45724620            sort ();
    45734621            for (;;) {
    … …  
    46084656                }
    46094657            }
    46104658        }
    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.
    46124660        iterator2 find2 (int rank, size_type i, size_type j, int direction = 1) {
    46134661            sort ();
    46144662            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}