Re: [boost] uBLAS: problem with sparse_vector

From: Joerg Walter (jhr.walter_at_[hidden])
Date: 2002-06-29 07:13:35

----- Original Message -----
From: "Marc Duflot" <m.duflot_at_[hidden]>
To: <boost_at_[hidden]>
Sent: Thursday, June 27, 2002 4:17 PM
Subject: [boost] uBLAS: problem with sparse_vector

> There is a bug in uBLAS that is illustrated with the following code:
>
> #include <iostream>
> #include <boost/numeric/ublas/vector_sp.h>
> #include <boost/numeric/ublas/io.h>
>
> int main () {
>     numerics::sparse_vector<double> v (3);
>     v(1) = v(0) = 42;
>     std::cout << v << std::endl;
> }
>
> The output of the program is (with g++ 2.95.2 and g++ 3.0.4 on Linux)
> [3](42,0,0)
>
> The element v(1) is 0 instead of 42. The bug comes from the fact that v(1)
is
> evaluated before v(0) = 42. The bug is not present with the Intel compiler
on
> the same platform.
>
> Note: there is the same problem with Loki::AssocVector.
>
> I think that this feature must be documented (or corrected if it is
possible)
> before uBLAS is accepted into boost.

I've just tried to fix this problem following Dave Abraham's hint (thanks,
Dave!). First of all I extended your sample to check std::map, map_array, ,
sparse_vector<..., std::map>, sparse_vector<..., map_array> and realized
this way, that the real problem was the semantic difference of std::map and
map_array.

Then I embedded the following nested proxy class in map_array:

        class proxy:
            public container_reference<map_array> {
        public:
            typedef typename map_array::index_type index_type;
            typedef typename map_array::data_value_type data_value_type;
            typedef const typename map_array::data_value_type
&data_const_reference;
            typedef typename map_array::data_value_type &data_reference;
            typedef std::pair<index_type, data_value_type> *pointer;

            // Construction and destruction
            proxy (map_array &a, pointer it):
                container_reference<map_array> (a), it_ (it), i_
(it->first), d_ (it->second) {
            }
            proxy (map_array &a, index_type i):
                container_reference<map_array> (a), i_ (i), d_ () {
                it_ = (*this) ().find (i_);
                if (it_ == (*this) ().end ())
                    it_ = (*this) ().insert ((*this) ().end (), value_type
(i_, d_));
                d_ = it_->second;
            }
            ~proxy () {
                if (it_->first != i_)
                    it_ = (*this) ().find (i_);
                it_->second = d_;
            }

            // Assignment
            template<class D>
            proxy &operator = (const D &d) {
                d_ = d;
                return *this;
            }
            proxy &operator = (const proxy &p) {
                d_ = p.d_;
                return *this;
            }
            template<class D>
            proxy &operator += (const D &d) {
                d_ += d;
                return *this;
            }
            proxy &operator += (const proxy &p) {
                d_ += p.d_;
                return *this;
            }
            template<class D>
            proxy &operator -= (const D &d) {
                d_ -= d;
                return *this;
            }
            proxy &operator -= (const proxy &p) {
                d_ -= p.d_;
                return *this;
            }
            template<class D>
            proxy &operator *= (const D &d) {
                d_ *= d;
                return *this;
            }
            proxy &operator *= (const proxy &p) {
                d_ *= p.d_;
                return *this;
            }
            template<class D>
            proxy &operator /= (const D &d) {
                d_ /= d;
                return *this;
            }
            proxy &operator /= (const proxy &p) {
                d_ /= p.d_;
                return *this;
            }

            // Conversion
            operator data_reference () {
                return d_;
            }

        private:
            pointer it_;
            index_type i_;
            data_value_type d_;
        };

and changed typedefs

        class proxy;
        typedef proxy data_reference;

and subscription operator

        // Element access
        NUMERICS_INLINE
        data_reference operator [] (index_type i) {
            // This fixes a [1] = a [0] = 1.
            // Thanks to Marc Duflot for spotting this.
            return data_reference (*this, i);
        }

With this changes, the test for map_array seemed to work (compiled with GCC
3.1), but: the interfaces of std::map and map_array are slightly different
now. So I
had to add the following map_traits

    template<class A>
    class map_traits {
        typedef void proxy;
    };
    template<class I, class T>
    class map_traits<std::map<I, T> > {
        typedef typename std::map<I, T>::data_type &proxy;
    };
    template<class I, class T>
    class map_traits<map_array<I, T> > {
        typedef typename map_array<I, T>::data_reference proxy;
    };

which use PTS (so the patch won't work on broken compilers) and then to
change in sparse_vector and sparse_matrix typedefs

        typedef typename map_traits<A>::proxy reference;

and mutable iterator dereference

            // Dereference
            NUMERICS_INLINE
            reference operator * () const {
                return reference ((*this) ().data (), it_);
            }

> Sorry if it is indeed documented somewhere or if it is a known bug.

May be, it would have been easier to only document that bug ;-)

Thanks for your feedback

Joerg