source: NonGTP/Boost/boost/numeric/ublas/vector_sparse.hpp @ 857

Revision 857, 73.6 KB checked in by igarcia, 18 years ago (diff)
Line 
1//
2//  Copyright (c) 2000-2002
3//  Joerg Walter, Mathias Koch
4//
5//  Permission to use, copy, modify, distribute and sell this software
6//  and its documentation for any purpose is hereby granted without fee,
7//  provided that the above copyright notice appear in all copies and
8//  that both that copyright notice and this permission notice appear
9//  in supporting documentation.  The authors make no representations
10//  about the suitability of this software for any purpose.
11//  It is provided "as is" without express or implied warranty.
12//
13//  The authors gratefully acknowledge the support of
14//  GeNeSys mbH & Co. KG in producing this work.
15//
16
17#ifndef _BOOST_UBLAS_VECTOR_SPARSE_
18#define _BOOST_UBLAS_VECTOR_SPARSE_
19
20#include <boost/numeric/ublas/storage_sparse.hpp>
21#include <boost/numeric/ublas/vector_expression.hpp>
22#include <boost/numeric/ublas/detail/vector_assign.hpp>
23
24// Iterators based on ideas of Jeremy Siek
25
26namespace boost { namespace numeric { namespace ublas {
27
28#ifdef BOOST_UBLAS_STRICT_VECTOR_SPARSE
29
30    template<class V>
31    class sparse_vector_element:
32       public container_reference<V> {
33    public:
34        typedef V vector_type;
35        typedef typename V::size_type size_type;
36        typedef typename V::value_type value_type;
37        typedef const value_type &const_reference;
38        typedef value_type *pointer;
39
40    private:
41        // Proxied element operations
42        void get_d () const {
43            pointer p = (*this) ().find_element (i_);
44            if (p)
45                d_ = *p;
46            else
47                d_ = value_type/*zero*/();
48        }
49
50        void set (const value_type &s) const {
51            pointer p = (*this) ().find_element (i_);
52            if (!p)
53                (*this) ().insert_element (i_, s);
54            else
55                *p = s;
56        }
57       
58    public:   
59        // Construction and destruction
60        sparse_vector_element (vector_type &v, size_type i):
61            container_reference<vector_type> (v), i_ (i) {
62        }
63        BOOST_UBLAS_INLINE
64        sparse_vector_element (const sparse_vector_element &p):
65            container_reference<vector_type> (p), i_ (p.i_) {}
66        BOOST_UBLAS_INLINE
67        ~sparse_vector_element () {
68        }
69
70        // Assignment
71        BOOST_UBLAS_INLINE
72        sparse_vector_element &operator = (const sparse_vector_element &p) {
73            // Overide the implict copy assignment
74            p.get_d ();
75            set (p.d_);
76            return *this;
77        }
78        template<class D>
79        BOOST_UBLAS_INLINE
80        sparse_vector_element &operator = (const D &d) {
81            set (d);
82            return *this;
83        }
84        template<class D>
85        BOOST_UBLAS_INLINE
86        sparse_vector_element &operator += (const D &d) {
87            get_d ();
88            d_ += d;
89            set (d_);
90            return *this;
91        }
92        template<class D>
93        BOOST_UBLAS_INLINE
94        sparse_vector_element &operator -= (const D &d) {
95            get_d ();
96            d_ -= d;
97            set (d_);
98            return *this;
99        }
100        template<class D>
101        BOOST_UBLAS_INLINE
102        sparse_vector_element &operator *= (const D &d) {
103            get_d ();
104            d_ *= d;
105            set (d_);
106            return *this;
107        }
108        template<class D>
109        BOOST_UBLAS_INLINE
110        sparse_vector_element &operator /= (const D &d) {
111            get_d ();
112            d_ /= d;
113            set (d_);
114            return *this;
115        }
116
117        // Comparison
118        template<class D>
119        BOOST_UBLAS_INLINE
120        bool operator == (const D &d) const {
121            get_d ();
122            return d_ == d;
123        }
124        template<class D>
125        BOOST_UBLAS_INLINE
126        bool operator != (const D &d) const {
127            get_d ();
128            return d_ != d;
129        }
130
131        // Conversion - weak link in proxy as d_ is not a perfect alias for the element
132        BOOST_UBLAS_INLINE
133        operator const_reference () const {
134            get_d ();
135            return d_;
136        }
137
138        // Conversion to reference - may be invalidated
139        BOOST_UBLAS_INLINE
140        value_type& ref () const {
141            pointer p = (*this) ().find_element (i_);
142            if (!p)
143                (*this) ().insert_element (i_, value_type/*zero*/());
144            return *p;
145        }
146
147    private:
148        size_type i_;
149        mutable value_type d_;
150    };
151
152    /*
153     * Generalise explicit reference access
154     */
155    namespace detail {
156        template <class R>
157        struct element_reference {
158            typedef R& reference;
159            static reference get_reference (reference r)
160            {
161                return r;
162            }
163        };
164        template <class V>
165        struct element_reference<sparse_vector_element<V> > {
166            typedef typename V::value_type& reference;
167            static reference get_reference (const sparse_vector_element<V>& sve)
168            {
169                return sve.ref ();
170            }
171        };
172    }
173    template <class VER>
174    typename detail::element_reference<VER>::reference ref (VER& ver) {
175        return detail::element_reference<VER>::get_reference (ver);
176    }
177    template <class VER>
178    typename detail::element_reference<VER>::reference ref (const VER& ver) {
179        return detail::element_reference<VER>::get_reference (ver);
180    }
181
182
183    template<class V>
184    struct type_traits<sparse_vector_element<V> > {
185        typedef typename V::value_type element_type;
186        typedef type_traits<sparse_vector_element<V> > self_type;
187        typedef typename type_traits<element_type>::value_type value_type;
188        typedef typename type_traits<element_type>::const_reference const_reference;
189        typedef sparse_vector_element<V> reference;
190        typedef typename type_traits<element_type>::real_type real_type;
191        typedef typename type_traits<element_type>::precision_type precision_type;
192
193        static const unsigned plus_complexity = type_traits<element_type>::plus_complexity;
194        static const unsigned multiplies_complexity = type_traits<element_type>::multiplies_complexity;
195
196        static
197        BOOST_UBLAS_INLINE
198        real_type real (const_reference t) {
199            return type_traits<element_type>::real (t);
200        }
201        static
202        BOOST_UBLAS_INLINE
203        real_type imag (const_reference t) {
204            return type_traits<element_type>::imag (t);
205        }
206        static
207        BOOST_UBLAS_INLINE
208        value_type conj (const_reference t) {
209            return type_traits<element_type>::conj (t);
210        }
211
212        static
213        BOOST_UBLAS_INLINE
214        real_type abs (const_reference t) {
215            return type_traits<element_type>::abs (t);
216        }
217        static
218        BOOST_UBLAS_INLINE
219        value_type sqrt (const_reference t) {
220            return type_traits<element_type>::sqrt (t);
221        }
222
223        static
224        BOOST_UBLAS_INLINE
225        real_type norm_1 (const_reference t) {
226            return type_traits<element_type>::norm_1 (t);
227        }
228        static
229        BOOST_UBLAS_INLINE
230        real_type norm_2 (const_reference t) {
231            return type_traits<element_type>::norm_2 (t);
232        }
233        static
234        BOOST_UBLAS_INLINE
235        real_type norm_inf (const_reference t) {
236            return type_traits<element_type>::norm_inf (t);
237        }
238
239        static
240        BOOST_UBLAS_INLINE
241        bool equals (const_reference t1, const_reference t2) {
242            return type_traits<element_type>::equals (t1, t2);
243        }
244    };
245
246    template<class V1, class T2>
247    struct promote_traits<sparse_vector_element<V1>, T2> {
248        typedef typename promote_traits<typename sparse_vector_element<V1>::value_type, T2>::promote_type promote_type;
249    };
250    template<class T1, class V2>
251    struct promote_traits<T1, sparse_vector_element<V2> > {
252        typedef typename promote_traits<T1, typename sparse_vector_element<V2>::value_type>::promote_type promote_type;
253    };
254    template<class V1, class V2>
255    struct promote_traits<sparse_vector_element<V1>, sparse_vector_element<V2> > {
256        typedef typename promote_traits<typename sparse_vector_element<V1>::value_type,
257                                        typename sparse_vector_element<V2>::value_type>::promote_type promote_type;
258    };
259
260#endif
261
262
263    // Index map based sparse vector class
264    template<class T, class A>
265    class mapped_vector:
266        public vector_container<mapped_vector<T, A> > {
267
268        typedef T &true_reference;
269        typedef T *pointer;
270        typedef const T *const_pointer;
271        typedef mapped_vector<T, A> self_type;
272    public:
273#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
274        using vector_container<self_type>::operator ();
275#endif
276        typedef typename A::size_type size_type;
277        typedef typename A::difference_type difference_type;
278        typedef T value_type;
279        typedef A array_type;
280        typedef const value_type &const_reference;
281#ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
282        typedef typename detail::map_traits<A,T>::reference reference;
283#else
284        typedef sparse_vector_element<self_type> reference;
285#endif
286        typedef const vector_reference<const self_type> const_closure_type;
287        typedef vector_reference<self_type> closure_type;
288        typedef self_type vector_temporary_type;
289        typedef sparse_tag storage_category;
290
291        // Construction and destruction
292        BOOST_UBLAS_INLINE
293        mapped_vector ():
294            vector_container<self_type> (),
295            size_ (0), data_ () {}
296        BOOST_UBLAS_INLINE
297        mapped_vector (size_type size, size_type non_zeros = 0):
298            vector_container<self_type> (),
299            size_ (size), data_ () {
300            detail::map_reserve (data(), restrict_capacity (non_zeros));
301        }
302        BOOST_UBLAS_INLINE
303        mapped_vector (const mapped_vector &v):
304            vector_container<self_type> (),
305            size_ (v.size_), data_ (v.data_) {}
306        template<class AE>
307        BOOST_UBLAS_INLINE
308        mapped_vector (const vector_expression<AE> &ae, size_type non_zeros = 0):
309            vector_container<self_type> (),
310            size_ (ae ().size ()), data_ () {
311            detail::map_reserve (data(), restrict_capacity (non_zeros));
312            vector_assign<scalar_assign> (*this, ae);
313        }
314
315        // Accessors
316        BOOST_UBLAS_INLINE
317        size_type size () const {
318            return size_;
319        }
320        BOOST_UBLAS_INLINE
321        size_type nnz_capacity () const {
322            return detail::map_capacity (data ());
323        }
324        BOOST_UBLAS_INLINE
325        size_type nnz () const {
326            return data (). size ();
327        }
328
329        // Storage accessors
330        BOOST_UBLAS_INLINE
331        const array_type &data () const {
332            return data_;
333        }
334        BOOST_UBLAS_INLINE
335        array_type &data () {
336            return data_;
337        }
338
339        // Resizing
340    private:
341        BOOST_UBLAS_INLINE
342        size_type restrict_capacity (size_type non_zeros) const {
343            non_zeros = (std::min) (non_zeros, size_);
344            return non_zeros;
345        }
346    public:
347        BOOST_UBLAS_INLINE
348        void resize (size_type size, bool preserve = true) {
349            size_ = size;
350            if (preserve) {
351                data ().erase (data ().lower_bound(size_), data ().end());
352            }
353            else {
354                data ().clear ();
355            }
356        }
357
358        // Reserving
359        BOOST_UBLAS_INLINE
360        void reserve (size_type non_zeros = 0, bool preserve = true) {
361            detail::map_reserve (data (), restrict_capacity (non_zeros));
362        }
363
364        // Element support
365        BOOST_UBLAS_INLINE
366        pointer find_element (size_type i) {
367            return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i));
368        }
369        BOOST_UBLAS_INLINE
370        const_pointer find_element (size_type i) const {
371            const_subiterator_type it (data ().find (i));
372            if (it == data ().end ())
373                return 0;
374            BOOST_UBLAS_CHECK ((*it).first == i, internal_logic ());   // broken map
375            return &(*it).second;
376        }
377
378        // Element access
379        BOOST_UBLAS_INLINE
380        const_reference operator () (size_type i) const {
381            BOOST_UBLAS_CHECK (i < size_, bad_index ());
382            const_subiterator_type it (data ().find (i));
383            if (it == data ().end ())
384                return zero_;
385            BOOST_UBLAS_CHECK ((*it).first == i, internal_logic ());   // broken map
386            return (*it).second;
387        }
388        BOOST_UBLAS_INLINE
389        true_reference ref (size_type i) {
390            BOOST_UBLAS_CHECK (i < size_, bad_index ());
391            std::pair<subiterator_type, bool> ii (data ().insert (typename array_type::value_type (i, value_type/*zero*/())));
392            BOOST_UBLAS_CHECK ((ii.first)->first == i, internal_logic ());   // broken map
393            return (ii.first)->second;
394        }
395        BOOST_UBLAS_INLINE
396        reference operator () (size_type i) {
397#ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
398            return ref (i);
399#else
400            BOOST_UBLAS_CHECK (i < size_, bad_index ());
401            return reference (*this, i);
402#endif
403        }
404
405        BOOST_UBLAS_INLINE
406        const_reference operator [] (size_type i) const {
407            return (*this) (i);
408        }
409        BOOST_UBLAS_INLINE
410        reference operator [] (size_type i) {
411            return (*this) (i);
412        }
413
414        // Element assignment
415        BOOST_UBLAS_INLINE
416        true_reference insert_element (size_type i, const_reference t) {
417            std::pair<subiterator_type, bool> ii = data ().insert (typename array_type::value_type (i, t));
418            BOOST_UBLAS_CHECK (ii.second, bad_index ());        // duplicate element
419            BOOST_UBLAS_CHECK ((ii.first)->first == i, internal_logic ());   // broken map
420            if (!ii.second)     // existing element
421                (ii.first)->second = t;
422            return (ii.first)->second;
423        }
424        BOOST_UBLAS_INLINE
425        void erase_element (size_type i) {
426            subiterator_type it = data ().find (i);
427            if (it == data ().end ())
428                return;
429            data ().erase (it);
430        }
431       
432        // Zeroing
433        BOOST_UBLAS_INLINE
434        void clear () {
435            data ().clear ();
436        }
437
438        // Assignment
439        BOOST_UBLAS_INLINE
440        mapped_vector &operator = (const mapped_vector &v) {
441            if (this != &v) {
442                size_ = v.size_;
443                data () = v.data ();
444            }
445            return *this;
446        }
447        template<class C>          // Container assignment without temporary
448        BOOST_UBLAS_INLINE
449        mapped_vector &operator = (const vector_container<C> &v) {
450            resize (v ().size (), false);
451            assign (v);
452            return *this;
453        }
454        BOOST_UBLAS_INLINE
455        mapped_vector &assign_temporary (mapped_vector &v) {
456            swap (v);
457            return *this;
458        }
459        template<class AE>
460        BOOST_UBLAS_INLINE
461        mapped_vector &operator = (const vector_expression<AE> &ae) {
462            self_type temporary (ae, detail::map_capacity (data()));
463            return assign_temporary (temporary);
464        }
465        template<class AE>
466        BOOST_UBLAS_INLINE
467        mapped_vector &assign (const vector_expression<AE> &ae) {
468            vector_assign<scalar_assign> (*this, ae);
469            return *this;
470        }
471
472        // Computed assignment
473        template<class AE>
474        BOOST_UBLAS_INLINE
475        mapped_vector &operator += (const vector_expression<AE> &ae) {
476            self_type temporary (*this + ae, detail::map_capacity (data()));
477            return assign_temporary (temporary);
478        }
479        template<class C>          // Container assignment without temporary
480        BOOST_UBLAS_INLINE
481        mapped_vector &operator += (const vector_container<C> &v) {
482            plus_assign (v);
483            return *this;
484        }
485        template<class AE>
486        BOOST_UBLAS_INLINE
487        mapped_vector &plus_assign (const vector_expression<AE> &ae) {
488            vector_assign<scalar_plus_assign> (*this, ae);
489            return *this;
490        }
491        template<class AE>
492        BOOST_UBLAS_INLINE
493        mapped_vector &operator -= (const vector_expression<AE> &ae) {
494            self_type temporary (*this - ae, detail::map_capacity (data()));
495            return assign_temporary (temporary);
496        }
497        template<class C>          // Container assignment without temporary
498        BOOST_UBLAS_INLINE
499        mapped_vector &operator -= (const vector_container<C> &v) {
500            minus_assign (v);
501            return *this;
502        }
503        template<class AE>
504        BOOST_UBLAS_INLINE
505        mapped_vector &minus_assign (const vector_expression<AE> &ae) {
506            vector_assign<scalar_minus_assign> (*this, ae);
507            return *this;
508        }
509        template<class AT>
510        BOOST_UBLAS_INLINE
511        mapped_vector &operator *= (const AT &at) {
512            vector_assign_scalar<scalar_multiplies_assign> (*this, at);
513            return *this;
514        }
515        template<class AT>
516        BOOST_UBLAS_INLINE
517        mapped_vector &operator /= (const AT &at) {
518            vector_assign_scalar<scalar_divides_assign> (*this, at);
519            return *this;
520        }
521
522        // Swapping
523        BOOST_UBLAS_INLINE
524        void swap (mapped_vector &v) {
525            if (this != &v) {
526                std::swap (size_, v.size_);
527                data ().swap (v.data ());
528            }
529        }
530        BOOST_UBLAS_INLINE
531        friend void swap (mapped_vector &v1, mapped_vector &v2) {
532            v1.swap (v2);
533        }
534
535        // Iterator types
536    private:
537        // Use storage iterator
538        typedef typename A::const_iterator const_subiterator_type;
539        typedef typename A::iterator subiterator_type;
540
541        BOOST_UBLAS_INLINE
542        true_reference at_element (size_type i) {
543            BOOST_UBLAS_CHECK (i < size_, bad_index ());
544            subiterator_type it (data ().find (i));
545            BOOST_UBLAS_CHECK (it != data ().end(), bad_index ());
546            BOOST_UBLAS_CHECK ((*it).first == i, internal_logic ());   // broken map
547            return it->second;
548        }
549
550    public:
551        class const_iterator;
552        class iterator;
553
554        // Element lookup
555        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
556        const_iterator find (size_type i) const {
557            return const_iterator (*this, data ().lower_bound (i));
558        }
559        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
560        iterator find (size_type i) {
561            return iterator (*this, data ().lower_bound (i));
562        }
563
564
565        class const_iterator:
566            public container_const_reference<mapped_vector>,
567            public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
568                                               const_iterator, value_type> {
569        public:
570            typedef typename mapped_vector::value_type value_type;
571            typedef typename mapped_vector::difference_type difference_type;
572            typedef typename mapped_vector::const_reference reference;
573            typedef const typename mapped_vector::pointer pointer;
574
575            // Construction and destruction
576            BOOST_UBLAS_INLINE
577            const_iterator ():
578                container_const_reference<self_type> (), it_ () {}
579            BOOST_UBLAS_INLINE
580            const_iterator (const self_type &v, const const_subiterator_type &it):
581                container_const_reference<self_type> (v), it_ (it) {}
582            BOOST_UBLAS_INLINE
583            const_iterator (const iterator &it):
584                container_const_reference<self_type> (it ()), it_ (it.it_) {}
585
586            // Arithmetic
587            BOOST_UBLAS_INLINE
588            const_iterator &operator ++ () {
589                ++ it_;
590                return *this;
591            }
592            BOOST_UBLAS_INLINE
593            const_iterator &operator -- () {
594                -- it_;
595                return *this;
596            }
597
598            // Dereference
599            BOOST_UBLAS_INLINE
600            const_reference operator * () const {
601                BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
602                return (*it_).second;
603            }
604
605            // Index
606            BOOST_UBLAS_INLINE
607            size_type index () const {
608                BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
609                BOOST_UBLAS_CHECK ((*it_).first < (*this) ().size (), bad_index ());
610                return (*it_).first;
611            }
612
613            // Assignment
614            BOOST_UBLAS_INLINE
615            const_iterator &operator = (const const_iterator &it) {
616                container_const_reference<self_type>::assign (&it ());
617                it_ = it.it_;
618                return *this;
619            }
620
621            // Comparison
622            BOOST_UBLAS_INLINE
623            bool operator == (const const_iterator &it) const {
624                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
625                return it_ == it.it_;
626            }
627
628        private:
629            const_subiterator_type it_;
630        };
631
632        BOOST_UBLAS_INLINE
633        const_iterator begin () const {
634            return const_iterator (*this, data ().begin ());
635        }
636        BOOST_UBLAS_INLINE
637        const_iterator end () const {
638            return const_iterator (*this, data ().end ());
639        }
640
641        class iterator:
642            public container_reference<mapped_vector>,
643            public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
644                                               iterator, value_type> {
645        public:
646            typedef typename mapped_vector::value_type value_type;
647            typedef typename mapped_vector::difference_type difference_type;
648            typedef typename mapped_vector::true_reference reference;
649            typedef typename mapped_vector::pointer pointer;
650
651            // Construction and destruction
652            BOOST_UBLAS_INLINE
653            iterator ():
654                container_reference<self_type> (), it_ () {}
655            BOOST_UBLAS_INLINE
656            iterator (self_type &v, const subiterator_type &it):
657                container_reference<self_type> (v), it_ (it) {}
658
659            // Arithmetic
660            BOOST_UBLAS_INLINE
661            iterator &operator ++ () {
662                ++ it_;
663                return *this;
664            }
665            BOOST_UBLAS_INLINE
666            iterator &operator -- () {
667                -- it_;
668                return *this;
669            }
670
671            // Dereference
672            BOOST_UBLAS_INLINE
673            reference operator * () const {
674                BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
675                return (*it_).second;
676            }
677
678            // Index
679            BOOST_UBLAS_INLINE
680            size_type index () const {
681                BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
682                BOOST_UBLAS_CHECK ((*it_).first < (*this) ().size (), bad_index ());
683                return (*it_).first;
684            }
685
686            // Assignment
687            BOOST_UBLAS_INLINE
688            iterator &operator = (const iterator &it) {
689                container_reference<self_type>::assign (&it ());
690                it_ = it.it_;
691                return *this;
692            }
693
694            // Comparison
695            BOOST_UBLAS_INLINE
696            bool operator == (const iterator &it) const {
697                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
698                return it_ == it.it_;
699            }
700
701        private:
702            subiterator_type it_;
703
704            friend class const_iterator;
705        };
706
707        BOOST_UBLAS_INLINE
708        iterator begin () {
709            return iterator (*this, data ().begin ());
710        }
711        BOOST_UBLAS_INLINE
712        iterator end () {
713            return iterator (*this, data ().end ());
714        }
715
716        // Reverse iterator
717        typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
718        typedef reverse_iterator_base<iterator> reverse_iterator;
719
720        BOOST_UBLAS_INLINE
721        const_reverse_iterator rbegin () const {
722            return const_reverse_iterator (end ());
723        }
724        BOOST_UBLAS_INLINE
725        const_reverse_iterator rend () const {
726            return const_reverse_iterator (begin ());
727        }
728        BOOST_UBLAS_INLINE
729        reverse_iterator rbegin () {
730            return reverse_iterator (end ());
731        }
732        BOOST_UBLAS_INLINE
733        reverse_iterator rend () {
734            return reverse_iterator (begin ());
735        }
736
737    private:
738        size_type size_;
739        array_type data_;
740        static const value_type zero_;
741    };
742
743    template<class T, class A>
744    const typename mapped_vector<T, A>::value_type mapped_vector<T, A>::zero_ = value_type/*zero*/();
745
746
747    // Compressed array based sparse vector class
748    // Thanks to Kresimir Fresl for extending this to cover different index bases.
749    template<class T, std::size_t IB, class IA, class TA>
750    class compressed_vector:
751        public vector_container<compressed_vector<T, IB, IA, TA> > {
752
753        typedef T &true_reference;
754        typedef T *pointer;
755        typedef const T *const_pointer;
756        typedef compressed_vector<T, IB, IA, TA> self_type;
757    public:
758#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
759        using vector_container<self_type>::operator ();
760#endif
761        // ISSUE require type consistency check
762        // is_convertable (IA::size_type, TA::size_type)
763        typedef typename IA::value_type size_type;
764        typedef typename IA::difference_type difference_type;
765        typedef T value_type;
766        typedef const T &const_reference;
767#ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
768        typedef T &reference;
769#else
770        typedef sparse_vector_element<self_type> reference;
771#endif
772        typedef IA index_array_type;
773        typedef TA value_array_type;
774        typedef const vector_reference<const self_type> const_closure_type;
775        typedef vector_reference<self_type> closure_type;
776        typedef self_type vector_temporary_type;
777        typedef sparse_tag storage_category;
778
779        // Construction and destruction
780        BOOST_UBLAS_INLINE
781        compressed_vector ():
782            vector_container<self_type> (),
783            size_ (0), capacity_ (restrict_capacity (0)), filled_ (0),
784            index_data_ (capacity_), value_data_ (capacity_) {
785            storage_invariants ();
786        }
787        explicit BOOST_UBLAS_INLINE
788        compressed_vector (size_type size, size_type non_zeros = 0):
789            vector_container<self_type> (),
790            size_ (size), capacity_ (restrict_capacity (non_zeros)), filled_ (0),
791            index_data_ (capacity_), value_data_ (capacity_) {
792        storage_invariants ();
793        }
794        BOOST_UBLAS_INLINE
795        compressed_vector (const compressed_vector &v):
796            vector_container<self_type> (),
797            size_ (v.size_), capacity_ (v.capacity_), filled_ (v.filled_),
798            index_data_ (v.index_data_), value_data_ (v.value_data_) {
799            storage_invariants ();
800        }
801        template<class AE>
802        BOOST_UBLAS_INLINE
803        compressed_vector (const vector_expression<AE> &ae, size_type non_zeros = 0):
804            vector_container<self_type> (),
805            size_ (ae ().size ()), capacity_ (restrict_capacity (non_zeros)), filled_ (0),
806            index_data_ (capacity_), value_data_ (capacity_) {
807            storage_invariants ();
808            vector_assign<scalar_assign> (*this, ae);
809        }
810
811        // Accessors
812        BOOST_UBLAS_INLINE
813        size_type size () const {
814            return size_;
815        }
816        BOOST_UBLAS_INLINE
817        size_type nnz_capacity () const {
818            return capacity_;
819        }
820        BOOST_UBLAS_INLINE
821        size_type nnz () const {
822            return filled_;
823        }
824
825        // Storage accessors
826        BOOST_UBLAS_INLINE
827        static size_type index_base () {
828            return IB;
829        }
830        BOOST_UBLAS_INLINE
831        typename index_array_type::size_type filled () const {
832            return filled_;
833        }
834        BOOST_UBLAS_INLINE
835        const index_array_type &index_data () const {
836            return index_data_;
837        }
838        BOOST_UBLAS_INLINE
839        const value_array_type &value_data () const {
840            return value_data_;
841        }
842        BOOST_UBLAS_INLINE
843        void set_filled (const typename index_array_type::size_type & filled) {
844            filled_ = filled;
845            storage_invariants ();
846        }
847        BOOST_UBLAS_INLINE
848        index_array_type &index_data () {
849            return index_data_;
850        }
851        BOOST_UBLAS_INLINE
852        value_array_type &value_data () {
853            return value_data_;
854        }
855
856        // Resizing
857    private:
858        BOOST_UBLAS_INLINE
859        size_type restrict_capacity (size_type non_zeros) const {
860            non_zeros = (std::max) (non_zeros, size_type (1));
861            non_zeros = (std::min) (non_zeros, size_);
862            return non_zeros;
863        }
864    public:
865        BOOST_UBLAS_INLINE
866        void resize (size_type size, bool preserve = true) {
867            // FIXME preserve unimplemented
868            BOOST_UBLAS_CHECK (!preserve, internal_logic ());
869            size_ = size;
870            capacity_ = restrict_capacity (capacity_);
871            index_data_. resize (capacity_);
872            value_data_. resize (capacity_);
873            filled_ = 0;
874            storage_invariants ();
875        }
876
877        // Reserving
878        BOOST_UBLAS_INLINE
879        void reserve (size_type non_zeros, bool preserve = true) {
880            capacity_ = restrict_capacity (non_zeros);
881            if (preserve) {
882                index_data_. resize (capacity_, size_type ());
883                value_data_. resize (capacity_, value_type ());
884                filled_ = (std::min) (capacity_, filled_);
885            }
886            else {
887                index_data_. resize (capacity_);
888                value_data_. resize (capacity_);
889                filled_ = 0;
890            }
891            storage_invariants ();
892        }
893
894        // Element support
895        BOOST_UBLAS_INLINE
896        pointer find_element (size_type i) {
897            return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i));
898        }
899        BOOST_UBLAS_INLINE
900        const_pointer find_element (size_type i) const {
901            const_subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
902            if (it == index_data_.begin () + filled_ || *it != k_based (i))
903                return 0;
904            return &value_data_ [it - index_data_.begin ()];
905        }
906
907        // Element access
908        BOOST_UBLAS_INLINE
909        const_reference operator () (size_type i) const {
910            BOOST_UBLAS_CHECK (i < size_, bad_index ());
911            const_subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
912            if (it == index_data_.begin () + filled_ || *it != k_based (i))
913                return zero_;
914            return value_data_ [it - index_data_.begin ()];
915        }
916        BOOST_UBLAS_INLINE
917        true_reference ref (size_type i) {
918            BOOST_UBLAS_CHECK (i < size_, bad_index ());
919            subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
920            if (it == index_data_.begin () + filled_ || *it != k_based (i))
921                return insert_element (i, value_type/*zero*/());
922            else
923                return value_data_ [it - index_data_.begin ()];
924        }
925        BOOST_UBLAS_INLINE
926        reference operator () (size_type i) {
927#ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
928            return ref (i) ;
929#else
930            BOOST_UBLAS_CHECK (i < size_, bad_index ());
931            return reference (*this, i);
932#endif
933        }
934
935        BOOST_UBLAS_INLINE
936        const_reference operator [] (size_type i) const {
937            return (*this) (i);
938        }
939        BOOST_UBLAS_INLINE
940        reference operator [] (size_type i) {
941            return (*this) (i);
942        }
943
944        // Element assignment
945        BOOST_UBLAS_INLINE
946        true_reference insert_element (size_type i, const_reference t) {
947            BOOST_UBLAS_CHECK (!find_element (i), bad_index ());        // duplicate element
948            if (filled_ >= capacity_)
949                reserve (2 * capacity_, true);
950            subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
951            // ISSUE max_capacity limit due to difference_type
952            typename std::iterator_traits<subiterator_type>::difference_type n = it - index_data_.begin ();
953            BOOST_UBLAS_CHECK (filled_ == 0 || filled_ == typename index_array_type::size_type (n) || *it != k_based (i), internal_logic ());   // duplicate found by lower_bound
954            ++ filled_;
955            it = index_data_.begin () + n;
956            std::copy_backward (it, index_data_.begin () + filled_ - 1, index_data_.begin () + filled_);
957            *it = k_based (i);
958            typename value_array_type::iterator itt (value_data_.begin () + n);
959            std::copy_backward (itt, value_data_.begin () + filled_ - 1, value_data_.begin () + filled_);
960            *itt = t;
961            storage_invariants ();
962            return *itt;
963        }
964        BOOST_UBLAS_INLINE
965        void erase_element (size_type i) {
966            subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
967            typename std::iterator_traits<subiterator_type>::difference_type  n = it - index_data_.begin ();
968            if (filled_ > typename index_array_type::size_type (n) && *it == k_based (i)) {
969                std::copy (it + 1, index_data_.begin () + filled_, it);
970                typename value_array_type::iterator itt (value_data_.begin () + n);
971                std::copy (itt + 1, value_data_.begin () + filled_, itt);
972                -- filled_;
973            }
974            storage_invariants ();
975        }
976       
977        // Zeroing
978        BOOST_UBLAS_INLINE
979        void clear () {
980            filled_ = 0;
981            storage_invariants ();
982        }
983
984        // Assignment
985        BOOST_UBLAS_INLINE
986        compressed_vector &operator = (const compressed_vector &v) {
987            if (this != &v) {
988                size_ = v.size_;
989                capacity_ = v.capacity_;
990                filled_ = v.filled_;
991                index_data_ = v.index_data_;
992                value_data_ = v.value_data_;
993            }
994            storage_invariants ();
995            return *this;
996        }
997        template<class C>          // Container assignment without temporary
998        BOOST_UBLAS_INLINE
999        compressed_vector &operator = (const vector_container<C> &v) {
1000            resize (v ().size (), false);
1001            assign (v);
1002            return *this;
1003        }
1004        BOOST_UBLAS_INLINE
1005        compressed_vector &assign_temporary (compressed_vector &v) {
1006            swap (v);
1007            return *this;
1008        }
1009        template<class AE>
1010        BOOST_UBLAS_INLINE
1011        compressed_vector &operator = (const vector_expression<AE> &ae) {
1012            self_type temporary (ae, capacity_);
1013            return assign_temporary (temporary);
1014        }
1015        template<class AE>
1016        BOOST_UBLAS_INLINE
1017        compressed_vector &assign (const vector_expression<AE> &ae) {
1018            vector_assign<scalar_assign> (*this, ae);
1019            return *this;
1020        }
1021
1022        // Computed assignment
1023        template<class AE>
1024        BOOST_UBLAS_INLINE
1025        compressed_vector &operator += (const vector_expression<AE> &ae) {
1026            self_type temporary (*this + ae, capacity_);
1027            return assign_temporary (temporary);
1028        }
1029        template<class C>          // Container assignment without temporary
1030        BOOST_UBLAS_INLINE
1031        compressed_vector &operator += (const vector_container<C> &v) {
1032            plus_assign (v);
1033            return *this;
1034        }
1035        template<class AE>
1036        BOOST_UBLAS_INLINE
1037        compressed_vector &plus_assign (const vector_expression<AE> &ae) {
1038            vector_assign<scalar_plus_assign> (*this, ae);
1039            return *this;
1040        }
1041        template<class AE>
1042        BOOST_UBLAS_INLINE
1043        compressed_vector &operator -= (const vector_expression<AE> &ae) {
1044            self_type temporary (*this - ae, capacity_);
1045            return assign_temporary (temporary);
1046        }
1047        template<class C>          // Container assignment without temporary
1048        BOOST_UBLAS_INLINE
1049        compressed_vector &operator -= (const vector_container<C> &v) {
1050            minus_assign (v);
1051            return *this;
1052        }
1053        template<class AE>
1054        BOOST_UBLAS_INLINE
1055        compressed_vector &minus_assign (const vector_expression<AE> &ae) {
1056            vector_assign<scalar_minus_assign> (*this, ae);
1057            return *this;
1058        }
1059        template<class AT>
1060        BOOST_UBLAS_INLINE
1061        compressed_vector &operator *= (const AT &at) {
1062            vector_assign_scalar<scalar_multiplies_assign> (*this, at);
1063            return *this;
1064        }
1065        template<class AT>
1066        BOOST_UBLAS_INLINE
1067        compressed_vector &operator /= (const AT &at) {
1068            vector_assign_scalar<scalar_divides_assign> (*this, at);
1069            return *this;
1070        }
1071
1072        // Swapping
1073        BOOST_UBLAS_INLINE
1074        void swap (compressed_vector &v) {
1075            if (this != &v) {
1076                std::swap (size_, v.size_);
1077                std::swap (capacity_, v.capacity_);
1078                std::swap (filled_, v.filled_);
1079                index_data_.swap (v.index_data_);
1080                value_data_.swap (v.value_data_);
1081            }
1082            storage_invariants ();
1083        }
1084        BOOST_UBLAS_INLINE
1085        friend void swap (compressed_vector &v1, compressed_vector &v2) {
1086            v1.swap (v2);
1087        }
1088
1089        // Back element insertion and erasure
1090        BOOST_UBLAS_INLINE
1091        void push_back (size_type i, const_reference t) {
1092            BOOST_UBLAS_CHECK (filled_ == 0 || index_data_ [filled_ - 1] < k_based (i), external_logic ());
1093            if (filled_ >= capacity_)
1094                reserve (2 * capacity_, true);
1095            BOOST_UBLAS_CHECK (filled_ < capacity_, internal_logic ());
1096            index_data_ [filled_] = k_based (i);
1097            value_data_ [filled_] = t;
1098            ++ filled_;
1099            storage_invariants ();
1100        }
1101        BOOST_UBLAS_INLINE
1102        void pop_back () {
1103            BOOST_UBLAS_CHECK (filled_ > 0, external_logic ());
1104            -- filled_;
1105            storage_invariants ();
1106        }
1107
1108        // Iterator types
1109    private:
1110        // Use index array iterator
1111        typedef typename IA::const_iterator const_subiterator_type;
1112        typedef typename IA::iterator subiterator_type;
1113
1114        BOOST_UBLAS_INLINE
1115        true_reference at_element (size_type i) {
1116            BOOST_UBLAS_CHECK (i < size_, bad_index ());
1117            subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1118            BOOST_UBLAS_CHECK (it != index_data_.begin () + filled_ && *it == k_based (i), bad_index ());
1119            return value_data_ [it - index_data_.begin ()];
1120        }
1121
1122    public:
1123        class const_iterator;
1124        class iterator;
1125
1126        // Element lookup
1127        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
1128        const_iterator find (size_type i) const {
1129            return const_iterator (*this, detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1130        }
1131        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
1132        iterator find (size_type i) {
1133            return iterator (*this, detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1134        }
1135
1136
1137        class const_iterator:
1138            public container_const_reference<compressed_vector>,
1139            public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
1140                                               const_iterator, value_type> {
1141        public:
1142            typedef typename compressed_vector::value_type value_type;
1143            typedef typename compressed_vector::difference_type difference_type;
1144            typedef typename compressed_vector::const_reference reference;
1145            typedef const typename compressed_vector::pointer pointer;
1146
1147            // Construction and destruction
1148            BOOST_UBLAS_INLINE
1149            const_iterator ():
1150                container_const_reference<self_type> (), it_ () {}
1151            BOOST_UBLAS_INLINE
1152            const_iterator (const self_type &v, const const_subiterator_type &it):
1153                container_const_reference<self_type> (v), it_ (it) {}
1154            BOOST_UBLAS_INLINE
1155            const_iterator (const iterator &it):
1156                container_const_reference<self_type> (it ()), it_ (it.it_) {}
1157
1158            // Arithmetic
1159            BOOST_UBLAS_INLINE
1160            const_iterator &operator ++ () {
1161                ++ it_;
1162                return *this;
1163            }
1164            BOOST_UBLAS_INLINE
1165            const_iterator &operator -- () {
1166                -- it_;
1167                return *this;
1168            }
1169
1170            // Dereference
1171            BOOST_UBLAS_INLINE
1172            const_reference operator * () const {
1173                BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
1174                return (*this) ().value_data_ [it_ - (*this) ().index_data_.begin ()];
1175            }
1176
1177            // Index
1178            BOOST_UBLAS_INLINE
1179            size_type index () const {
1180                BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
1181                BOOST_UBLAS_CHECK ((*this) ().zero_based (*it_) < (*this) ().size (), bad_index ());
1182                return (*this) ().zero_based (*it_);
1183            }
1184
1185            // Assignment
1186            BOOST_UBLAS_INLINE
1187            const_iterator &operator = (const const_iterator &it) {
1188                container_const_reference<self_type>::assign (&it ());
1189                it_ = it.it_;
1190                return *this;
1191            }
1192
1193            // Comparison
1194            BOOST_UBLAS_INLINE
1195            bool operator == (const const_iterator &it) const {
1196                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1197                return it_ == it.it_;
1198            }
1199
1200        private:
1201            const_subiterator_type it_;
1202        };
1203
1204        BOOST_UBLAS_INLINE
1205        const_iterator begin () const {
1206            return find (0);
1207        }
1208        BOOST_UBLAS_INLINE
1209        const_iterator end () const {
1210            return find (size_);
1211        }
1212
1213        class iterator:
1214            public container_reference<compressed_vector>,
1215            public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
1216                                               iterator, value_type> {
1217        public:
1218            typedef typename compressed_vector::value_type value_type;
1219            typedef typename compressed_vector::difference_type difference_type;
1220            typedef typename compressed_vector::true_reference reference;
1221            typedef typename compressed_vector::pointer pointer;
1222
1223            // Construction and destruction
1224            BOOST_UBLAS_INLINE
1225            iterator ():
1226                container_reference<self_type> (), it_ () {}
1227            BOOST_UBLAS_INLINE
1228            iterator (self_type &v, const subiterator_type &it):
1229                container_reference<self_type> (v), it_ (it) {}
1230
1231            // Arithmetic
1232            BOOST_UBLAS_INLINE
1233            iterator &operator ++ () {
1234                ++ it_;
1235                return *this;
1236            }
1237            BOOST_UBLAS_INLINE
1238            iterator &operator -- () {
1239                -- it_;
1240                return *this;
1241            }
1242
1243            // Dereference
1244            BOOST_UBLAS_INLINE
1245            reference operator * () const {
1246                BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
1247                return (*this) ().value_data_ [it_ - (*this) ().index_data_.begin ()];
1248            }
1249
1250            // Index
1251            BOOST_UBLAS_INLINE
1252            size_type index () const {
1253                BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
1254                BOOST_UBLAS_CHECK ((*this) ().zero_based (*it_) < (*this) ().size (), bad_index ());
1255                return (*this) ().zero_based (*it_);
1256            }
1257
1258            // Assignment
1259            BOOST_UBLAS_INLINE
1260            iterator &operator = (const iterator &it) {
1261                container_reference<self_type>::assign (&it ());
1262                it_ = it.it_;
1263                return *this;
1264            }
1265
1266            // Comparison
1267            BOOST_UBLAS_INLINE
1268            bool operator == (const iterator &it) const {
1269                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1270                return it_ == it.it_;
1271            }
1272
1273        private:
1274            subiterator_type it_;
1275
1276            friend class const_iterator;
1277        };
1278
1279        BOOST_UBLAS_INLINE
1280        iterator begin () {
1281            return find (0);
1282        }
1283        BOOST_UBLAS_INLINE
1284        iterator end () {
1285            return find (size_);
1286        }
1287
1288        // Reverse iterator
1289        typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
1290        typedef reverse_iterator_base<iterator> reverse_iterator;
1291
1292        BOOST_UBLAS_INLINE
1293        const_reverse_iterator rbegin () const {
1294            return const_reverse_iterator (end ());
1295        }
1296        BOOST_UBLAS_INLINE
1297        const_reverse_iterator rend () const {
1298            return const_reverse_iterator (begin ());
1299        }
1300        BOOST_UBLAS_INLINE
1301        reverse_iterator rbegin () {
1302            return reverse_iterator (end ());
1303        }
1304        BOOST_UBLAS_INLINE
1305        reverse_iterator rend () {
1306            return reverse_iterator (begin ());
1307        }
1308
1309    private:
1310        void storage_invariants () const
1311        {
1312            BOOST_UBLAS_CHECK (capacity_ == index_data_.size (), internal_logic ());
1313            BOOST_UBLAS_CHECK (capacity_ == value_data_.size (), internal_logic ());
1314            BOOST_UBLAS_CHECK (filled_ <= capacity_, internal_logic ());
1315        }
1316
1317        size_type size_;
1318        typename index_array_type::size_type capacity_;
1319        typename index_array_type::size_type filled_;
1320        index_array_type index_data_;
1321        value_array_type value_data_;
1322        static const value_type zero_;
1323
1324        BOOST_UBLAS_INLINE
1325        static size_type zero_based (size_type k_based_index) {
1326            return k_based_index - IB;
1327        }
1328        BOOST_UBLAS_INLINE
1329        static size_type k_based (size_type zero_based_index) {
1330            return zero_based_index + IB;
1331        }
1332
1333        friend class iterator;
1334        friend class const_iterator;
1335    };
1336
1337    template<class T, std::size_t IB, class IA, class TA>
1338    const typename compressed_vector<T, IB, IA, TA>::value_type compressed_vector<T, IB, IA, TA>::zero_ = value_type/*zero*/();
1339
1340
1341    // Coordimate array based sparse vector class
1342    // Thanks to Kresimir Fresl for extending this to cover different index bases.
1343    template<class T, std::size_t IB, class IA, class TA>
1344    class coordinate_vector:
1345        public vector_container<coordinate_vector<T, IB, IA, TA> > {
1346
1347        typedef T &true_reference;
1348        typedef T *pointer;
1349        typedef const T *const_pointer;
1350        typedef coordinate_vector<T, IB, IA, TA> self_type;
1351    public:
1352#ifdef BOOST_UBLAS_ENABLE_PROXY_SHORTCUTS
1353        using vector_container<self_type>::operator ();
1354#endif
1355        // ISSUE require type consistency check
1356        // is_convertable (IA::size_type, TA::size_type)
1357        typedef typename IA::value_type size_type;
1358        typedef typename IA::difference_type difference_type;
1359        typedef T value_type;
1360        typedef const T &const_reference;
1361#ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
1362        typedef T &reference;
1363#else
1364        typedef sparse_vector_element<self_type> reference;
1365#endif
1366        typedef IA index_array_type;
1367        typedef TA value_array_type;
1368        typedef const vector_reference<const self_type> const_closure_type;
1369        typedef vector_reference<self_type> closure_type;
1370        typedef self_type vector_temporary_type;
1371        typedef sparse_tag storage_category;
1372
1373        // Construction and destruction
1374        BOOST_UBLAS_INLINE
1375        coordinate_vector ():
1376            vector_container<self_type> (),
1377            size_ (0), capacity_ (restrict_capacity (0)),
1378            filled_ (0), sorted_filled_ (filled_), sorted_ (true),
1379            index_data_ (capacity_), value_data_ (capacity_) {
1380            storage_invariants ();
1381        }
1382        explicit BOOST_UBLAS_INLINE
1383        coordinate_vector (size_type size, size_type non_zeros = 0):
1384            vector_container<self_type> (),
1385            size_ (size), capacity_ (restrict_capacity (non_zeros)),
1386            filled_ (0), sorted_filled_ (filled_), sorted_ (true),
1387            index_data_ (capacity_), value_data_ (capacity_) {
1388            storage_invariants ();
1389        }
1390        BOOST_UBLAS_INLINE
1391        coordinate_vector (const coordinate_vector &v):
1392            vector_container<self_type> (),
1393            size_ (v.size_), capacity_ (v.capacity_),
1394            filled_ (v.filled_), sorted_filled_ (v.sorted_filled_), sorted_ (v.sorted_),
1395            index_data_ (v.index_data_), value_data_ (v.value_data_) {
1396            storage_invariants ();
1397        }
1398        template<class AE>
1399        BOOST_UBLAS_INLINE
1400        coordinate_vector (const vector_expression<AE> &ae, size_type non_zeros = 0):
1401            vector_container<self_type> (),
1402            size_ (ae ().size ()), capacity_ (restrict_capacity (non_zeros)),
1403            filled_ (0), sorted_filled_ (filled_), sorted_ (true),
1404            index_data_ (capacity_), value_data_ (capacity_) {
1405            storage_invariants ();
1406            vector_assign<scalar_assign> (*this, ae);
1407        }
1408
1409        // Accessors
1410        BOOST_UBLAS_INLINE
1411        size_type size () const {
1412            return size_;
1413        }
1414        BOOST_UBLAS_INLINE
1415        size_type nnz_capacity () const {
1416            return capacity_;
1417        }
1418        BOOST_UBLAS_INLINE
1419        size_type nnz () const {
1420            return filled_;
1421        }
1422
1423        // Storage accessors
1424        BOOST_UBLAS_INLINE
1425        static size_type index_base () {
1426            return IB;
1427        }
1428        BOOST_UBLAS_INLINE
1429        typename index_array_type::size_type filled () const {
1430            return filled_;
1431        }
1432        BOOST_UBLAS_INLINE
1433        const index_array_type &index_data () const {
1434            return index_data_;
1435        }
1436        BOOST_UBLAS_INLINE
1437        const value_array_type &value_data () const {
1438            return value_data_;
1439        }
1440        BOOST_UBLAS_INLINE
1441        void set_filled (const typename index_array_type::size_type &sorted, const typename index_array_type::size_type &filled) {
1442            sorted_filled_ = sorted;
1443            filled_ = filled;
1444            storage_invariants ();
1445            return filled_;
1446        }
1447        BOOST_UBLAS_INLINE
1448        index_array_type &index_data () {
1449            return index_data_;
1450        }
1451        BOOST_UBLAS_INLINE
1452        value_array_type &value_data () {
1453            return value_data_;
1454        }
1455
1456        // Resizing
1457    private:
1458        BOOST_UBLAS_INLINE
1459        size_type restrict_capacity (size_type non_zeros) const {
1460             // minimum non_zeros
1461             non_zeros = (std::max) (non_zeros, size_type (1));
1462             // ISSUE no maximum as coordinate may contain inserted duplicates
1463             return non_zeros;
1464        }
1465    public:
1466        BOOST_UBLAS_INLINE
1467        void resize (size_type size, bool preserve = true) {
1468            if (preserve)
1469                sort ();    // remove duplicate elements.
1470            capacity_ = restrict_capacity (capacity_);
1471            if (preserve) {
1472                index_data_. resize (capacity_, size_type ());
1473                value_data_. resize (capacity_, value_type ());
1474                filled_ = (std::min) (capacity_, filled_);
1475            }
1476            else {
1477                index_data_. resize (capacity_);
1478                value_data_. resize (capacity_);
1479                filled_ = 0;
1480            }
1481            sorted_filled_ = filled_;
1482            size_ = size;
1483            storage_invariants ();
1484        }
1485        // Reserving
1486        BOOST_UBLAS_INLINE
1487        void reserve (size_type non_zeros, bool preserve = true) {
1488            if (preserve)
1489                sort ();    // remove duplicate elements.
1490            capacity_ = restrict_capacity (non_zeros);
1491            if (preserve) {
1492                index_data_. resize (capacity_, size_type ());
1493                value_data_. resize (capacity_, value_type ());
1494                filled_ = (std::min) (capacity_, filled_);
1495                }
1496            else {
1497                index_data_. resize (capacity_);
1498                value_data_. resize (capacity_);
1499                filled_ = 0;
1500            }
1501            sorted_filled_ = filled_;
1502            storage_invariants ();
1503        }
1504
1505        // Element support
1506        BOOST_UBLAS_INLINE
1507        pointer find_element (size_type i) {
1508            return const_cast<pointer> (const_cast<const self_type&>(*this).find_element (i));
1509        }
1510        BOOST_UBLAS_INLINE
1511        const_pointer find_element (size_type i) const {
1512            sort ();
1513            const_subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1514            if (it == index_data_.begin () + filled_ || *it != k_based (i))
1515                return 0;
1516            return &value_data_ [it - index_data_.begin ()];
1517        }
1518
1519        // Element access
1520        BOOST_UBLAS_INLINE
1521        const_reference operator () (size_type i) const {
1522            BOOST_UBLAS_CHECK (i < size_, bad_index ());
1523            sort ();
1524            const_subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1525            if (it == index_data_.begin () + filled_ || *it != k_based (i))
1526                return zero_;
1527            return value_data_ [it - index_data_.begin ()];
1528        }
1529        BOOST_UBLAS_INLINE
1530        true_reference ref (size_type i) {
1531            BOOST_UBLAS_CHECK (i < size_, bad_index ());
1532            sort ();
1533            subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1534            if (it == index_data_.begin () + filled_ || *it != k_based (i))
1535                return insert_element (i, value_type/*zero*/());
1536            else
1537                return value_data_ [it - index_data_.begin ()];
1538        }
1539        BOOST_UBLAS_INLINE
1540        reference operator () (size_type i) {
1541#ifndef BOOST_UBLAS_STRICT_VECTOR_SPARSE
1542            return ref (i);
1543#else
1544            BOOST_UBLAS_CHECK (i < size_, bad_index ());
1545            return reference (*this, i);
1546#endif
1547        }
1548
1549        BOOST_UBLAS_INLINE
1550        const_reference operator [] (size_type i) const {
1551            return (*this) (i);
1552        }
1553        BOOST_UBLAS_INLINE
1554        reference operator [] (size_type i) {
1555            return (*this) (i);
1556        }
1557
1558        // Element assignment
1559        BOOST_UBLAS_INLINE
1560        void append_element (size_type i, const_reference t) {
1561            if (filled_ >= capacity_)
1562                reserve (2 * filled_, true);
1563            BOOST_UBLAS_CHECK (filled_ < capacity_, internal_logic ());
1564            index_data_ [filled_] = k_based (i);
1565            value_data_ [filled_] = t;
1566            ++ filled_;
1567            sorted_ = false;
1568            storage_invariants ();
1569        }
1570        BOOST_UBLAS_INLINE
1571        true_reference insert_element (size_type i, const_reference t) {
1572            BOOST_UBLAS_CHECK (!find_element (i), bad_index ());        // duplicate element
1573            append_element (i, t);
1574            return value_data_ [filled_ - 1];
1575        }
1576        BOOST_UBLAS_INLINE
1577        void erase_element (size_type i) {
1578            sort ();
1579            subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1580            typename std::iterator_traits<subiterator_type>::difference_type n = it - index_data_.begin ();
1581            if (filled_ > typename index_array_type::size_type (n) && *it == k_based (i)) {
1582                std::copy (it + 1, index_data_.begin () + filled_, it);
1583                typename value_array_type::iterator itt (value_data_.begin () + n);
1584                std::copy (itt + 1, value_data_.begin () + filled_, itt);
1585                -- filled_;
1586                sorted_filled_ = filled_;
1587            }
1588            storage_invariants ();
1589        }
1590       
1591        // Zeroing
1592        BOOST_UBLAS_INLINE
1593        void clear () {
1594            filled_ = 0;
1595            sorted_filled_ = filled_;
1596            sorted_ = true;
1597            storage_invariants ();
1598        }
1599       
1600        // Assignment
1601        BOOST_UBLAS_INLINE
1602        coordinate_vector &operator = (const coordinate_vector &v) {
1603            if (this != &v) {
1604                size_ = v.size_;
1605                capacity_ = v.capacity_;
1606                filled_ = v.filled_;
1607                sorted_filled_ = v.sorted_filled_;
1608                sorted_ = v.sorted_;
1609                index_data_ = v.index_data_;
1610                value_data_ = v.value_data_;
1611            }
1612            storage_invariants ();
1613            return *this;
1614        }
1615        template<class C>          // Container assignment without temporary
1616        BOOST_UBLAS_INLINE
1617        coordinate_vector &operator = (const vector_container<C> &v) {
1618            resize (v ().size (), false);
1619            assign (v);
1620            return *this;
1621        }
1622        BOOST_UBLAS_INLINE
1623        coordinate_vector &assign_temporary (coordinate_vector &v) {
1624            swap (v);
1625            return *this;
1626        }
1627        template<class AE>
1628        BOOST_UBLAS_INLINE
1629        coordinate_vector &operator = (const vector_expression<AE> &ae) {
1630            self_type temporary (ae, capacity_);
1631            return assign_temporary (temporary);
1632        }
1633        template<class AE>
1634        BOOST_UBLAS_INLINE
1635        coordinate_vector &assign (const vector_expression<AE> &ae) {
1636            vector_assign<scalar_assign> (*this, ae);
1637            return *this;
1638        }
1639
1640        // Computed assignment
1641        template<class AE>
1642        BOOST_UBLAS_INLINE
1643        coordinate_vector &operator += (const vector_expression<AE> &ae) {
1644            self_type temporary (*this + ae, capacity_);
1645            return assign_temporary (temporary);
1646        }
1647        template<class C>          // Container assignment without temporary
1648        BOOST_UBLAS_INLINE
1649        coordinate_vector &operator += (const vector_container<C> &v) {
1650            plus_assign (v);
1651            return *this;
1652        }
1653        template<class AE>
1654        BOOST_UBLAS_INLINE
1655        coordinate_vector &plus_assign (const vector_expression<AE> &ae) {
1656            vector_assign<scalar_plus_assign> (*this, ae);
1657            return *this;
1658        }
1659        template<class AE>
1660        BOOST_UBLAS_INLINE
1661        coordinate_vector &operator -= (const vector_expression<AE> &ae) {
1662            self_type temporary (*this - ae, capacity_);
1663            return assign_temporary (temporary);
1664        }
1665        template<class C>          // Container assignment without temporary
1666        BOOST_UBLAS_INLINE
1667        coordinate_vector &operator -= (const vector_container<C> &v) {
1668            minus_assign (v);
1669            return *this;
1670        }
1671        template<class AE>
1672        BOOST_UBLAS_INLINE
1673        coordinate_vector &minus_assign (const vector_expression<AE> &ae) {
1674            vector_assign<scalar_minus_assign> (*this, ae);
1675            return *this;
1676        }
1677        template<class AT>
1678        BOOST_UBLAS_INLINE
1679        coordinate_vector &operator *= (const AT &at) {
1680            vector_assign_scalar<scalar_multiplies_assign> (*this, at);
1681            return *this;
1682        }
1683        template<class AT>
1684        BOOST_UBLAS_INLINE
1685        coordinate_vector &operator /= (const AT &at) {
1686            vector_assign_scalar<scalar_divides_assign> (*this, at);
1687            return *this;
1688        }
1689
1690        // Swapping
1691        BOOST_UBLAS_INLINE
1692        void swap (coordinate_vector &v) {
1693            if (this != &v) {
1694                std::swap (size_, v.size_);
1695                std::swap (capacity_, v.capacity_);
1696                std::swap (filled_, v.filled_);
1697                std::swap (sorted_filled_, v.sorted_filled_);
1698                std::swap (sorted_, v.sorted_);
1699                index_data_.swap (v.index_data_);
1700                value_data_.swap (v.value_data_);
1701            }
1702            storage_invariants ();
1703        }
1704        BOOST_UBLAS_INLINE
1705        friend void swap (coordinate_vector &v1, coordinate_vector &v2) {
1706            v1.swap (v2);
1707        }
1708
1709        // Sorting and summation of duplicates
1710        BOOST_UBLAS_INLINE
1711        void sort () const {
1712            if (! sorted_ && filled_ > 0) {
1713                typedef index_pair_array<index_array_type, value_array_type> array_pair;
1714                array_pair ipa (filled_, index_data_, value_data_);
1715                const typename array_pair::iterator iunsorted = ipa.begin () + sorted_filled_;
1716                // sort new elements and merge
1717                std::sort (iunsorted, ipa.end ());
1718                std::inplace_merge (ipa.begin (), iunsorted, ipa.end ());
1719               
1720                // sum duplicates with += and remove
1721                size_type filled = 0;
1722                for (size_type i = 1; i < filled_; ++ i) {
1723                    if (index_data_ [filled] != index_data_ [i]) {
1724                        ++ filled;
1725                        if (filled != i) {
1726                            index_data_ [filled] = index_data_ [i];
1727                            value_data_ [filled] = value_data_ [i];
1728                        }
1729                    } else {
1730                        value_data_ [filled] += value_data_ [i];
1731                    }
1732                }
1733                filled_ = filled + 1;
1734                sorted_filled_ = filled_;
1735                sorted_ = true;
1736                storage_invariants ();
1737            }
1738        }
1739
1740        // Back element insertion and erasure
1741        BOOST_UBLAS_INLINE
1742        void push_back (size_type i, const_reference t) {
1743            // must maintain sort order
1744            BOOST_UBLAS_CHECK (sorted_ && (filled_ == 0 || index_data_ [filled_ - 1] < k_based (i)), external_logic ());
1745            if (filled_ >= capacity_)
1746                reserve (2 * filled_, true);
1747            BOOST_UBLAS_CHECK (filled_ < capacity_, internal_logic ());
1748            index_data_ [filled_] = k_based (i);
1749            value_data_ [filled_] = t;
1750            ++ filled_;
1751            sorted_filled_ = filled_;
1752            storage_invariants ();
1753        }
1754        BOOST_UBLAS_INLINE
1755        void pop_back () {
1756            // ISSUE invariants could be simpilfied if sorted required as precondition
1757            BOOST_UBLAS_CHECK (filled_ > 0, external_logic ());
1758            -- filled_;
1759            sorted_filled_ = (std::min) (sorted_filled_, filled_);
1760            sorted_ = sorted_filled_ = filled_;
1761            storage_invariants ();
1762        }
1763
1764        // Iterator types
1765    private:
1766        // Use index array iterator
1767        typedef typename IA::const_iterator const_subiterator_type;
1768        typedef typename IA::iterator subiterator_type;
1769
1770        BOOST_UBLAS_INLINE
1771        true_reference at_element (size_type i) {
1772            BOOST_UBLAS_CHECK (i < size_, bad_index ());
1773            sort ();
1774            subiterator_type it (detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1775            BOOST_UBLAS_CHECK (it != index_data_.begin () + filled_ && *it == k_based (i), bad_index ());
1776            return value_data_ [it - index_data_.begin ()];
1777        }
1778
1779    public:
1780        class const_iterator;
1781        class iterator;
1782
1783        // Element lookup
1784        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
1785        const_iterator find (size_type i) const {
1786            sort ();
1787            return const_iterator (*this, detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1788        }
1789        // BOOST_UBLAS_INLINE This function seems to be big. So we do not let the compiler inline it.   
1790        iterator find (size_type i) {
1791            sort ();
1792            return iterator (*this, detail::lower_bound (index_data_.begin (), index_data_.begin () + filled_, k_based (i), std::less<size_type> ()));
1793        }
1794
1795
1796        class const_iterator:
1797            public container_const_reference<coordinate_vector>,
1798            public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
1799                                               const_iterator, value_type> {
1800        public:
1801            typedef typename coordinate_vector::value_type value_type;
1802            typedef typename coordinate_vector::difference_type difference_type;
1803            typedef typename coordinate_vector::const_reference reference;
1804            typedef const typename coordinate_vector::pointer pointer;
1805
1806            // Construction and destruction
1807            BOOST_UBLAS_INLINE
1808            const_iterator ():
1809                container_const_reference<self_type> (), it_ () {}
1810            BOOST_UBLAS_INLINE
1811            const_iterator (const self_type &v, const const_subiterator_type &it):
1812                container_const_reference<self_type> (v), it_ (it) {}
1813            BOOST_UBLAS_INLINE
1814            const_iterator (const iterator &it):
1815                container_const_reference<self_type> (it ()), it_ (it.it_) {}
1816
1817            // Arithmetic
1818            BOOST_UBLAS_INLINE
1819            const_iterator &operator ++ () {
1820                ++ it_;
1821                return *this;
1822            }
1823            BOOST_UBLAS_INLINE
1824            const_iterator &operator -- () {
1825                -- it_;
1826                return *this;
1827            }
1828
1829            // Dereference
1830            BOOST_UBLAS_INLINE
1831            const_reference operator * () const {
1832                BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
1833                return (*this) ().value_data_ [it_ - (*this) ().index_data_.begin ()];
1834            }
1835
1836            // Index
1837            BOOST_UBLAS_INLINE
1838            size_type index () const {
1839                BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
1840                BOOST_UBLAS_CHECK ((*this) ().zero_based (*it_) < (*this) ().size (), bad_index ());
1841                return (*this) ().zero_based (*it_);
1842            }
1843
1844            // Assignment
1845            BOOST_UBLAS_INLINE
1846            const_iterator &operator = (const const_iterator &it) {
1847                container_const_reference<self_type>::assign (&it ());
1848                it_ = it.it_;
1849                return *this;
1850            }
1851
1852            // Comparison
1853            BOOST_UBLAS_INLINE
1854            bool operator == (const const_iterator &it) const {
1855                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1856                return it_ == it.it_;
1857            }
1858
1859        private:
1860            const_subiterator_type it_;
1861        };
1862
1863        BOOST_UBLAS_INLINE
1864        const_iterator begin () const {
1865            return find (0);
1866        }
1867        BOOST_UBLAS_INLINE
1868        const_iterator end () const {
1869            return find (size_);
1870        }
1871
1872        class iterator:
1873            public container_reference<coordinate_vector>,
1874            public bidirectional_iterator_base<sparse_bidirectional_iterator_tag,
1875                                               iterator, value_type> {
1876        public:
1877            typedef typename coordinate_vector::value_type value_type;
1878            typedef typename coordinate_vector::difference_type difference_type;
1879            typedef typename coordinate_vector::true_reference reference;
1880            typedef typename coordinate_vector::pointer pointer;
1881
1882            // Construction and destruction
1883            BOOST_UBLAS_INLINE
1884            iterator ():
1885                container_reference<self_type> (), it_ () {}
1886            BOOST_UBLAS_INLINE
1887            iterator (self_type &v, const subiterator_type &it):
1888                container_reference<self_type> (v), it_ (it) {}
1889
1890            // Arithmetic
1891            BOOST_UBLAS_INLINE
1892            iterator &operator ++ () {
1893                ++ it_;
1894                return *this;
1895            }
1896            BOOST_UBLAS_INLINE
1897            iterator &operator -- () {
1898                -- it_;
1899                return *this;
1900            }
1901
1902            // Dereference
1903            BOOST_UBLAS_INLINE
1904            reference operator * () const {
1905                BOOST_UBLAS_CHECK (index () < (*this) ().size (), bad_index ());
1906                return (*this) ().value_data_ [it_ - (*this) ().index_data_.begin ()];
1907            }
1908
1909            // Index
1910            BOOST_UBLAS_INLINE
1911            size_type index () const {
1912                BOOST_UBLAS_CHECK (*this != (*this) ().end (), bad_index ());
1913                BOOST_UBLAS_CHECK ((*this) ().zero_based (*it_) < (*this) ().size (), bad_index ());
1914                return (*this) ().zero_based (*it_);
1915            }
1916
1917            // Assignment
1918            BOOST_UBLAS_INLINE
1919            iterator &operator = (const iterator &it) {
1920                container_reference<self_type>::assign (&it ());
1921                it_ = it.it_;
1922                return *this;
1923            }
1924
1925            // Comparison
1926            BOOST_UBLAS_INLINE
1927            bool operator == (const iterator &it) const {
1928                BOOST_UBLAS_CHECK (&(*this) () == &it (), external_logic ());
1929                return it_ == it.it_;
1930            }
1931
1932        private:
1933            subiterator_type it_;
1934
1935            friend class const_iterator;
1936        };
1937
1938        BOOST_UBLAS_INLINE
1939        iterator begin () {
1940            return find (0);
1941        }
1942        BOOST_UBLAS_INLINE
1943        iterator end () {
1944            return find (size_);
1945        }
1946
1947        // Reverse iterator
1948        typedef reverse_iterator_base<const_iterator> const_reverse_iterator;
1949        typedef reverse_iterator_base<iterator> reverse_iterator;
1950
1951        BOOST_UBLAS_INLINE
1952        const_reverse_iterator rbegin () const {
1953            return const_reverse_iterator (end ());
1954        }
1955        BOOST_UBLAS_INLINE
1956        const_reverse_iterator rend () const {
1957            return const_reverse_iterator (begin ());
1958        }
1959        BOOST_UBLAS_INLINE
1960        reverse_iterator rbegin () {
1961            return reverse_iterator (end ());
1962        }
1963        BOOST_UBLAS_INLINE
1964        reverse_iterator rend () {
1965            return reverse_iterator (begin ());
1966        }
1967
1968    private:
1969        void storage_invariants () const
1970        {
1971            BOOST_UBLAS_CHECK (capacity_ == index_data_.size (), internal_logic ());
1972            BOOST_UBLAS_CHECK (capacity_ == value_data_.size (), internal_logic ());
1973            BOOST_UBLAS_CHECK (filled_ <= capacity_, internal_logic ());
1974            BOOST_UBLAS_CHECK (sorted_filled_ <= filled_, internal_logic ());
1975            BOOST_UBLAS_CHECK (sorted_ == (sorted_filled_ == filled_), internal_logic ());
1976        }
1977
1978        size_type size_;
1979        size_type capacity_;
1980        mutable typename index_array_type::size_type filled_;
1981        mutable typename index_array_type::size_type sorted_filled_;
1982        mutable bool sorted_;
1983        mutable index_array_type index_data_;
1984        mutable value_array_type value_data_;
1985        static const value_type zero_;
1986
1987        BOOST_UBLAS_INLINE
1988        static size_type zero_based (size_type k_based_index) {
1989            return k_based_index - IB;
1990        }
1991        BOOST_UBLAS_INLINE
1992        static size_type k_based (size_type zero_based_index) {
1993            return zero_based_index + IB;
1994        }
1995
1996        friend class iterator;
1997        friend class const_iterator;
1998    };
1999
2000    template<class T, std::size_t IB, class IA, class TA>
2001    const typename coordinate_vector<T, IB, IA, TA>::value_type coordinate_vector<T, IB, IA, TA>::zero_ = value_type/*zero*/();
2002
2003}}}
2004
2005#endif
Note: See TracBrowser for help on using the repository browser.