[boost] Re: ublas & custom containers? (2)

From: Kresimir Fresl (fresl_at_[hidden])
Date: 2002-03-26 06:04:50

Joerg Walter wrote:

> [...] If a sparse matrix result is assigned to a dense matrix,
 > the assignment operation should be dispatched to
 > the sparse matrix logic: clear the dense matrix and insert
> the rhs elements at the corresponding positions.

I am afraid that I found a bug in the implementation of this
assignment:

This works:

    // default array_type: unbounded_array<>
    typedef numerics::vector<double> vct_t;
    typedef numerics::sparse_vector<double> sp_vct_t;

    vct_t v (8);
    sp_vct_t sv (8, 3);

    sv[1] = 1.1; sv[3] = 3.1; sv[4] = 4.1;

    v = sv;
    cout << v << endl;
    // output is, of course: [8](0,1.1,0,3.1,4.1,0,0,0)

This works, too:

    typedef numerics::vector< double,
                 numerics::array_adaptor<double> > vct_t;
    typedef numerics::sparse_vector<double> sp_vct_t;

    std::vector<double> storage (8);   // ;o)
    numerics::array_adaptor<double> aa (8, &storage[0]);
    vct_t v (8, aa);

    vct_t v (8);
    sp_vct_t sv (8, 3);
    sv[1] = 1.1; sv[3] = 3.1; sv[4] = 4.1;

    v = sv; // same output

But this doesn't:

    typedef numerics::vector< double,
                                              std::vector<double> > vct_t;
    typedef numerics::sparse_vector<double> sp_vct_t;

    vct_t v (8);
    sp_vct_t sv (8, 3);
    sv[1] = 1.1; sv[3] = 3.1; sv[4] = 4.1;

    v = sv;

It causes ``Segmentation fault'' (at least on Linux
with g++ 3.0.4).  And IMHO with good reason:

`numerics::vector_assign<>::operator()' with `sparse_tag'
uses `array_type::clear()' to clear `numerics::vector<>'
before assignement, that is, before the call to
`array_type::insert()'. But the behaviour and semantics of
`numerics::unbounded_array<>::clear()' and
`numerics::array_adaptor<>::clear()' is very different
from the behaviour of `std::vector<>::clear()'.

First two simply assign `0' to all elements. But
`std::vector<>::clear()' really clears the vector,
ie. it *removes* all elements, and the postcondition
is, according to standard, `size() == 0' which
means that the vector is empty. And after that
the insertation at particular position (except at the
very beginning) fails because this position simply
doesn't exist.

Either the lines

      if (v.data().size() == 0)
         v.data().resize (e().size());

should be added after `v.clear()' in
`vector_assign<>::operator() (vct, vct_expr, sparse_tag)',
or

      std::fill (v.begin(), v.end(), value_type());

should be used instead of `v.clear ();'.

Another possibilities are to specialize `vector_assign<>'
for `std::vector<>' or to use some traits class.

Sincerely,

fres

PS. BTW, arguments of the constructor of `sparse_vector<>'
are `size' and `non_zero', that is vector size and number of
non-zeros. Documentation says only: ``Allocates a
sparse_vector that holds at most size elements.'' What about
`non_zero'? Can it hold more than `non_zero' non-zeros?
Less? Or must one call `resize (size, non_zero2)' to change
the number of non-zeros?

PPS. It would be useful to add `num_non_zeros()' function
to `sparse_vector<>' which returns the number of non-zero
elements.