Boost.uBlas 1.49
Linear Algebra in C++: matrices, vectors and numeric algorithms

storage_sparse.hpp

Go to the documentation of this file.
00001 //
00002 //  Copyright (c) 2000-2002
00003 //  Joerg Walter, Mathias Koch
00004 //
00005 //  Distributed under the Boost Software License, Version 1.0. (See
00006 //  accompanying file LICENSE_1_0.txt or copy at
00007 //  http://www.boost.org/LICENSE_1_0.txt)
00008 //
00009 //  The authors gratefully acknowledge the support of
00010 //  GeNeSys mbH & Co. KG in producing this work.
00011 //
00012 
00013 #ifndef _BOOST_UBLAS_STORAGE_SPARSE_
00014 #define _BOOST_UBLAS_STORAGE_SPARSE_
00015 
00016 #include <map>
00017 #include <boost/serialization/collection_size_type.hpp>
00018 #include <boost/serialization/nvp.hpp>
00019 #include <boost/serialization/array.hpp>
00020 #include <boost/serialization/map.hpp>
00021 #include <boost/serialization/base_object.hpp>
00022 
00023 #include <boost/numeric/ublas/storage.hpp>
00024 
00025 
00026 namespace boost { namespace numeric { namespace ublas {
00027 
00028     namespace detail {
00029 
00030         template<class I, class T, class C>
00031         BOOST_UBLAS_INLINE
00032         I lower_bound (const I &begin, const I &end, const T &t, C compare) {
00033             // t <= *begin <=> ! (*begin < t)
00034             if (begin == end || ! compare (*begin, t))
00035                 return begin;
00036             if (compare (*(end - 1), t))
00037                 return end;
00038             return std::lower_bound (begin, end, t, compare);
00039         }
00040         template<class I, class T, class C>
00041         BOOST_UBLAS_INLINE
00042         I upper_bound (const I &begin, const I &end, const T &t, C compare) {
00043             if (begin == end || compare (t, *begin))
00044                 return begin;
00045             // (*end - 1) <= t <=> ! (t < *end)
00046             if (! compare (t, *(end - 1)))
00047                 return end;
00048             return std::upper_bound (begin, end, t, compare);
00049         }
00050 
00051         template<class P>
00052         struct less_pair {
00053             BOOST_UBLAS_INLINE
00054             bool operator () (const P &p1, const P &p2) {
00055                 return p1.first < p2.first;
00056             }
00057         };
00058         template<class T>
00059         struct less_triple {
00060             BOOST_UBLAS_INLINE
00061             bool operator () (const T &t1, const T &t2) {
00062                 return t1.first.first < t2.first.first ||
00063                        (t1.first.first == t2.first.first && t1.first.second < t2.first.second);
00064             }
00065         };
00066 
00067     }
00068 
00069 #ifdef BOOST_UBLAS_STRICT_MAP_ARRAY
00070     template<class A>
00071     class sparse_storage_element:
00072        public container_reference<A> {
00073     public:
00074         typedef A array_type;
00075         typedef typename A::key_type index_type;
00076         typedef typename A::mapped_type data_value_type;
00077         // typedef const data_value_type &data_const_reference;
00078         typedef typename type_traits<data_value_type>::const_reference data_const_reference;
00079         typedef data_value_type &data_reference;
00080         typedef typename A::value_type value_type;
00081         typedef value_type *pointer;
00082 
00083         // Construction and destruction
00084         BOOST_UBLAS_INLINE
00085         sparse_storage_element (array_type &a, pointer it):
00086             container_reference<array_type> (a), it_ (it), i_ (it->first), d_ (it->second), dirty_ (false) {}
00087         BOOST_UBLAS_INLINE
00088         sparse_storage_element (array_type &a, index_type i):
00089             container_reference<array_type> (a), it_ (), i_ (i), d_ (), dirty_ (false) {
00090             pointer it = (*this) ().find (i_);
00091             if (it == (*this) ().end ())
00092                 it = (*this) ().insert ((*this) ().end (), value_type (i_, d_));
00093             d_ = it->second;
00094         }
00095         BOOST_UBLAS_INLINE
00096         ~sparse_storage_element () {
00097             if (dirty_) {
00098                 if (! it_)
00099                     it_ = (*this) ().find (i_);
00100                 BOOST_UBLAS_CHECK (it_ != (*this) ().end (), internal_logic ());
00101                 it_->second = d_;
00102             }
00103         }
00104 
00105         // Element access - only if data_const_reference is defined
00106         BOOST_UBLAS_INLINE
00107         typename data_value_type::data_const_reference
00108         operator [] (index_type i) const {
00109             return d_ [i];
00110         }
00111 
00112         // Assignment
00113         BOOST_UBLAS_INLINE
00114         sparse_storage_element &operator = (const sparse_storage_element &p) {
00115             // Overide the implict copy assignment
00116             d_ = p.d_;
00117             dirty_ = true;
00118             return *this;
00119         }
00120         template<class D>
00121         BOOST_UBLAS_INLINE
00122         sparse_storage_element &operator = (const D &d) {
00123             d_ = d;
00124             dirty_ = true;
00125             return *this;
00126         }
00127         template<class D>
00128         BOOST_UBLAS_INLINE
00129         sparse_storage_element &operator += (const D &d) {
00130             d_ += d;
00131             dirty_ = true;
00132             return *this;
00133         }
00134         template<class D>
00135         BOOST_UBLAS_INLINE
00136         sparse_storage_element &operator -= (const D &d) {
00137             d_ -= d;
00138             dirty_ = true;
00139             return *this;
00140         }
00141         template<class D>
00142         BOOST_UBLAS_INLINE
00143         sparse_storage_element &operator *= (const D &d) {
00144             d_ *= d;
00145             dirty_ = true;
00146             return *this;
00147         }
00148         template<class D>
00149         BOOST_UBLAS_INLINE
00150         sparse_storage_element &operator /= (const D &d) {
00151             d_ /= d;
00152             dirty_ = true;
00153             return *this;
00154         }
00155 
00156         // Comparison
00157         template<class D>
00158         BOOST_UBLAS_INLINE
00159         bool operator == (const D &d) const {
00160             return d_ == d;
00161         }
00162         template<class D>
00163         BOOST_UBLAS_INLINE
00164         bool operator != (const D &d) const {
00165             return d_ != d;
00166         }
00167 
00168         // Conversion
00169         BOOST_UBLAS_INLINE
00170         operator data_const_reference () const {
00171             return d_;
00172         }
00173 
00174         // Swapping
00175         BOOST_UBLAS_INLINE
00176         void swap (sparse_storage_element p) {
00177             if (this != &p) {
00178                 dirty_ = true;
00179                 p.dirty_ = true;
00180                 std::swap (d_, p.d_);
00181             }
00182         }
00183         BOOST_UBLAS_INLINE
00184         friend void swap (sparse_storage_element p1, sparse_storage_element p2) {
00185             p1.swap (p2);
00186         }
00187 
00188     private:
00189         pointer it_;
00190         index_type i_;
00191         data_value_type d_;
00192         bool dirty_;
00193     };
00194 #endif
00195 
00196 
00197     // Default map type is simply forwarded to std::map
00198     // FIXME should use ALLOC for map but std::allocator of std::pair<const I, T> and std::pair<I,T> fail to compile
00199     template<class I, class T, class ALLOC>
00200     class map_std : public std::map<I, T /*, ALLOC */> {
00201     public:
00202          // Serialization
00203         template<class Archive>
00204         void serialize(Archive & ar, const unsigned int /* file_version */){
00205             ar & serialization::make_nvp("base", boost::serialization::base_object< std::map<I, T /*, ALLOC */> >(*this));
00206         }
00207     };
00208 
00209     
00210 
00211 
00212     // Map array
00213     //  Implementation requires pair<I, T> allocator definition (without const)
00214     template<class I, class T, class ALLOC>
00215     class map_array {
00216     public:
00217         typedef ALLOC allocator_type;
00218         typedef typename ALLOC::size_type size_type;
00219         typedef typename ALLOC::difference_type difference_type;
00220         typedef std::pair<I,T> value_type;
00221         typedef I key_type;
00222         typedef T mapped_type;
00223         typedef const value_type &const_reference;
00224         typedef value_type &reference;
00225         typedef const value_type *const_pointer;
00226         typedef value_type *pointer;
00227         // Iterators simply are pointers.
00228         typedef const_pointer const_iterator;
00229         typedef pointer iterator;
00230 
00231         typedef const T &data_const_reference;
00232 #ifndef BOOST_UBLAS_STRICT_MAP_ARRAY
00233         typedef T &data_reference;
00234 #else
00235         typedef sparse_storage_element<map_array> data_reference;
00236 #endif
00237 
00238         // Construction and destruction
00239         BOOST_UBLAS_INLINE
00240         map_array (const ALLOC &a = ALLOC()):
00241             alloc_(a), capacity_ (0), size_ (0) {
00242                 data_ = 0;
00243         }
00244         BOOST_UBLAS_INLINE
00245         map_array (const map_array &c):
00246             alloc_ (c.alloc_), capacity_ (c.size_), size_ (c.size_) {
00247             if (capacity_) {
00248                 data_ = alloc_.allocate (capacity_);
00249                 std::uninitialized_copy (data_, data_ + capacity_, c.data_);
00250                 // capacity != size_ requires uninitialized_fill (size_ to capacity_)
00251             }
00252             else
00253                 data_ = 0;
00254         }
00255         BOOST_UBLAS_INLINE
00256         ~map_array () {
00257             if (capacity_) {
00258                 std::for_each (data_, data_ + capacity_, static_destroy);
00259                 alloc_.deallocate (data_, capacity_);
00260             }
00261         }
00262 
00263     private:
00264         // Resizing - implicitly exposses uninitialized (but default constructed) mapped_type
00265         BOOST_UBLAS_INLINE
00266         void resize (size_type size) {
00267             BOOST_UBLAS_CHECK (size_ <= capacity_, internal_logic ());
00268             if (size > capacity_) {
00269                 const size_type capacity = size << 1;
00270                 BOOST_UBLAS_CHECK (capacity, internal_logic ());
00271                 pointer data = alloc_.allocate (capacity);
00272                 std::uninitialized_copy (data_, data_ + (std::min) (size, size_), data);
00273                 std::uninitialized_fill (data + (std::min) (size, size_), data + capacity, value_type ());
00274 
00275                 if (capacity_) {
00276                     std::for_each (data_, data_ + capacity_, static_destroy);
00277                     alloc_.deallocate (data_, capacity_);
00278                 }
00279                 capacity_ = capacity;
00280                 data_ = data;
00281             }
00282             size_ = size;
00283             BOOST_UBLAS_CHECK (size_ <= capacity_, internal_logic ());
00284         }
00285     public:
00286 
00287         // Reserving
00288         BOOST_UBLAS_INLINE
00289         void reserve (size_type capacity) {
00290             BOOST_UBLAS_CHECK (size_ <= capacity_, internal_logic ());
00291             // Reduce capacity_ if size_ allows
00292             BOOST_UBLAS_CHECK (capacity >= size_, bad_size ());
00293             pointer data;
00294             if (capacity) {
00295                 data = alloc_.allocate (capacity);
00296                 std::uninitialized_copy (data_, data_ + size_, data);
00297                 std::uninitialized_fill (data + size_, data + capacity, value_type ());
00298             }
00299             else
00300                 data = 0;
00301                 
00302             if (capacity_) {
00303                 std::for_each (data_, data_ + capacity_, static_destroy);
00304                 alloc_.deallocate (data_, capacity_);
00305             }
00306             capacity_ = capacity;
00307             data_ = data;
00308             BOOST_UBLAS_CHECK (size_ <= capacity_, internal_logic ());
00309         }
00310 
00311         // Random Access Container
00312         BOOST_UBLAS_INLINE
00313         size_type size () const {
00314             return size_;
00315         }
00316         BOOST_UBLAS_INLINE
00317         size_type capacity () const {
00318             return capacity_;
00319         }
00320         BOOST_UBLAS_INLINE
00321         size_type max_size () const {
00322             return 0; //TODO
00323         }
00324        
00325         BOOST_UBLAS_INLINE
00326         bool empty () const {
00327             return size_ == 0;
00328         }
00329             
00330         // Element access
00331         BOOST_UBLAS_INLINE
00332         data_reference operator [] (key_type i) {
00333 #ifndef BOOST_UBLAS_STRICT_MAP_ARRAY
00334             pointer it = find (i);
00335             if (it == end ())
00336                 it = insert (end (), value_type (i, mapped_type (0)));
00337             BOOST_UBLAS_CHECK (it != end (), internal_logic ());
00338             return it->second;
00339 #else
00340             return data_reference (*this, i);
00341 #endif
00342         }
00343 
00344         // Assignment
00345         BOOST_UBLAS_INLINE
00346         map_array &operator = (const map_array &a) {
00347             if (this != &a) {
00348                 resize (a.size_);
00349                 std::copy (a.data_, a.data_ + a.size_, data_);
00350             }
00351             return *this;
00352         }
00353         BOOST_UBLAS_INLINE
00354         map_array &assign_temporary (map_array &a) {
00355             swap (a);
00356             return *this;
00357         }
00358 
00359         // Swapping
00360         BOOST_UBLAS_INLINE
00361         void swap (map_array &a) {
00362             if (this != &a) {
00363                 std::swap (capacity_, a.capacity_);
00364                 std::swap (data_, a.data_);
00365                 std::swap (size_, a.size_);
00366             }
00367         }
00368         BOOST_UBLAS_INLINE
00369         friend void swap (map_array &a1, map_array &a2) {
00370             a1.swap (a2);
00371         }
00372 
00373         // Element insertion and deletion
00374         
00375         // From Back Insertion Sequence concept
00376         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
00377         iterator push_back (iterator it, const value_type &p) {
00378             if (size () == 0 || (it = end () - 1)->first < p.first) {
00379                 resize (size () + 1);
00380                 *(it = end () - 1) = p;
00381                 return it;
00382             }
00383             external_logic ().raise ();
00384             return it;
00385         }
00386         // Form Unique Associative Container concept
00387         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
00388         std::pair<iterator,bool> insert (const value_type &p) {
00389             iterator it = detail::lower_bound (begin (), end (), p, detail::less_pair<value_type> ());
00390             if (it != end () && it->first == p.first)
00391                 return std::make_pair (it, false);
00392             difference_type n = it - begin ();
00393             resize (size () + 1);
00394             it = begin () + n;    // allow for invalidation
00395             std::copy_backward (it, end () - 1, end ());
00396             *it = p;
00397             return std::make_pair (it, true);
00398         }
00399         // Form Sorted Associative Container concept
00400         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
00401         iterator insert (iterator hint, const value_type &p) {
00402             return insert (p).first;
00403         }
00404         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
00405         void erase (iterator it) {
00406             BOOST_UBLAS_CHECK (begin () <= it && it < end (), bad_index ());
00407             std::copy (it + 1, end (), it);
00408             resize (size () - 1);
00409         }
00410         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
00411         void erase (iterator it1, iterator it2) {
00412             if (it1 == it2) return /* nothing to erase */;
00413             BOOST_UBLAS_CHECK (begin () <= it1 && it1 < it2 && it2 <= end (), bad_index ());
00414             std::copy (it2, end (), it1);
00415             resize (size () - (it2 - it1));
00416         }
00417         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
00418         void clear () {
00419             resize (0);
00420         }
00421 
00422         // Element lookup
00423         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
00424         const_iterator find (key_type i) const {
00425             const_iterator it (detail::lower_bound (begin (), end (), value_type (i, mapped_type (0)), detail::less_pair<value_type> ()));
00426             if (it == end () || it->first != i)
00427                 it = end ();
00428             return it;
00429         }
00430         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
00431         iterator find (key_type i) {
00432             iterator it (detail::lower_bound (begin (), end (), value_type (i, mapped_type (0)), detail::less_pair<value_type> ()));
00433             if (it == end () || it->first != i)
00434                 it = end ();
00435             return it;
00436         }
00437         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
00438         const_iterator lower_bound (key_type i) const {
00439             return detail::lower_bound (begin (), end (), value_type (i, mapped_type (0)), detail::less_pair<value_type> ());
00440         }
00441         // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.    
00442         iterator lower_bound (key_type i) {
00443             return detail::lower_bound (begin (), end (), value_type (i, mapped_type (0)), detail::less_pair<value_type> ());
00444         }
00445 
00446         BOOST_UBLAS_INLINE
00447         const_iterator begin () const {
00448             return data_;
00449         }
00450         BOOST_UBLAS_INLINE
00451         const_iterator end () const {
00452             return data_ + size_;
00453         }
00454 
00455         BOOST_UBLAS_INLINE
00456         iterator begin () {
00457             return data_;
00458         }
00459         BOOST_UBLAS_INLINE
00460         iterator end () {
00461             return data_ + size_;
00462         }
00463 
00464         // Reverse iterators
00465         typedef std::reverse_iterator<const_iterator> const_reverse_iterator;
00466         typedef std::reverse_iterator<iterator> reverse_iterator;
00467 
00468         BOOST_UBLAS_INLINE
00469         const_reverse_iterator rbegin () const {
00470             return const_reverse_iterator (end ());
00471         }
00472         BOOST_UBLAS_INLINE
00473         const_reverse_iterator rend () const {
00474             return const_reverse_iterator (begin ());
00475         }
00476         BOOST_UBLAS_INLINE
00477         reverse_iterator rbegin () {
00478             return reverse_iterator (end ());
00479         }
00480         BOOST_UBLAS_INLINE
00481         reverse_iterator rend () {
00482             return reverse_iterator (begin ());
00483         }
00484 
00485         // Allocator
00486         allocator_type get_allocator () {
00487             return alloc_;
00488         }
00489 
00490          // Serialization
00491         template<class Archive>
00492         void serialize(Archive & ar, const unsigned int /* file_version */){
00493             serialization::collection_size_type s (size_);
00494             ar & serialization::make_nvp("size",s);
00495             if (Archive::is_loading::value) {
00496                 resize(s);
00497             }
00498             ar & serialization::make_array(data_, s);
00499         }
00500 
00501     private:
00502         // Provide destroy as a non member function
00503         BOOST_UBLAS_INLINE
00504         static void static_destroy (reference p) {
00505             (&p) -> ~value_type ();
00506         }
00507         ALLOC alloc_;
00508         size_type capacity_;
00509         pointer data_;
00510         size_type size_;
00511     };
00512 
00513 
00514     namespace detail {
00515         template<class A, class T>
00516         struct map_traits {
00517             typedef typename A::mapped_type &reference;
00518         };
00519         template<class I, class T, class ALLOC>
00520         struct map_traits<map_array<I, T, ALLOC>, T > {
00521             typedef typename map_array<I, T, ALLOC>::data_reference reference;
00522         };
00523 
00524         // reserve helpers for map_array and generic maps
00525         // ISSUE should be in map_traits but want to use on all compilers
00526 
00527         template<class M>
00528         BOOST_UBLAS_INLINE
00529         void map_reserve (M &/* m */, typename M::size_type /* capacity */) {
00530         }
00531         template<class I, class T, class ALLOC>
00532         BOOST_UBLAS_INLINE
00533         void map_reserve (map_array<I, T, ALLOC> &m, typename map_array<I, T, ALLOC>::size_type capacity) {
00534             m.reserve (capacity);
00535         }
00536 
00537         template<class M>
00538         struct map_capacity_traits {
00539             typedef typename M::size_type type ;
00540             type operator() ( M const& m ) const {
00541                return m.size ();
00542             }
00543         } ;
00544 
00545         template<class I, class T, class ALLOC>
00546         struct map_capacity_traits< map_array<I, T, ALLOC> > {
00547             typedef typename map_array<I, T, ALLOC>::size_type type ;
00548             type operator() ( map_array<I, T, ALLOC> const& m ) const {
00549                return m.capacity ();
00550             }
00551         } ;
00552 
00553         template<class M>
00554         BOOST_UBLAS_INLINE
00555         typename map_capacity_traits<M>::type map_capacity (M const& m) {
00556             return map_capacity_traits<M>() ( m );
00557         }
00558     }
00559 
00560 }}}
00561 
00562 #endif