Re: [boost] ublas: iterator problem in sparse matrix

From: Joerg Walter (jhr.walter_at_[hidden])
Date: 2002-06-24 14:31:29

----- Original Message -----
From: "Benedikt Weber" <weber_at_[hidden]>
Newsgroups: gmane.comp.lib.boost.devel
To: <boost_at_[hidden]>
Sent: Sunday, June 23, 2002 8:11 PM
Subject: [boost] ublas: iterator problem in sparse matrix

> I was playing around with ublas today to better understand how it works.
One
> point that got my attention is the 2-dimensional iterator (officially
> "iterator1" and "iterator2"). I was assuming that the iterators of a
sparse
> matrix iterate over the non-zero elements,

Correct.

> but when trying it out, I found
> the following:
>
> numerics::sparse_matrix<double> A(3,3,1);
>  A(0,0)= 1;
>
>  for (numerics::sparse_matrix<double>::iterator1 i= A.begin1(); i !=
> A.end1(); ++i)
>   for (numerics::sparse_matrix<double>::iterator2 j= i.begin(); j !=
> i.end(); ++j)
>     std::cout << "(" << i.index1() << "," << j.index2() << ") ";
>  std::cout << std::endl;
>
>  This prints out:
>  (0,0) (1,0) (1,1) (1,2) (2,0) (2,1) (2,2)
>
> Is this intention and what would the reason, or is it a bug?

I confirm, that this is a bug.

I slightly extended your test case:

    numerics::sparse_matrix<double> A(3,3,1);
    A.clear();
    A(0,0)= 1;
    for (numerics::sparse_matrix<double>::iterator1 i= A.begin1(); i !=
A.end1(); ++i)
        for (numerics::sparse_matrix<double>::iterator2 j= i.begin(); j !=
i.end(); ++j)
                std::cout << "(" << i.index1() << "," << j.index2() << ") ";
    std::cout << std::endl;
    A.clear();
    A(2,2)= 1;
    for (numerics::sparse_matrix<double>::iterator1 i= A.begin1(); i !=
A.end1(); ++i)
        for (numerics::sparse_matrix<double>::iterator2 j= i.begin(); j !=
i.end(); ++j)
                std::cout << "(" << i.index1() << "," << j.index2() << ") ";
    std::cout << std::endl;
    numerics::compressed_matrix<double> B(3,3,1);
    B.clear();
    B(0,0)= 1;
    for (numerics::compressed_matrix<double>::iterator1 i= B.begin1(); i !=
B.end1(); ++i)
        for (numerics::compressed_matrix<double>::iterator2 j= i.begin(); j
!= i.end(); ++j)
                std::cout << "(" << i.index1() << "," << j.index2() << ") ";
    std::cout << std::endl;
    B.clear();
    B(2,2)= 1;
    for (numerics::compressed_matrix<double>::iterator1 i= B.begin1(); i !=
B.end1(); ++i)
        for (numerics::compressed_matrix<double>::iterator2 j= i.begin(); j
!= i.end(); ++j)
                std::cout << "(" << i.index1() << "," << j.index2() << ") ";
    std::cout << std::endl;

All other cases gave the expected result ;-(

> The consequence is that matrix multiplication takes very long for sparse
> matrices. The following code takes forever, even though there is only one
> multiplication involved.
>
> // large matrices with just one entry
>  numerics::sparse_matrix<double> B(1000,1000,1);
>  numerics::sparse_matrix<double> C(1000,1000,1);
>  B(0,0)= 1;
>  C(0,0)= 1;
> numerics::sparse_matrix<double> D= numerics::prod(B,C);

Not surprising ;-(

> I certainly hope this is an easy to fix implementation bug or a wrong
> interpretation on my side and not a design problem. (sparse_vector works a
> expected and so is the sparse_matrix/sparce_vector multiplication)

The following bug fix to sparse_matrix<> seems to work for me:

        // Element lookup
        const_iterator1 find1 (int rank, size_type i, size_type j) const {
            const_iterator_type it (data ().lower_bound
(functor_type::element (i, size1_, j, size2_)));
            const_iterator_type it_end (data ().end ());
            size_type index2 = size_type (-1);
            while (it != it_end) {
                index2 = functor_type::index2 ((*it).first, size1_, size2_);
                if ((rank == 0 && index2 >= j) ||
                    (rank == 1 && index2 == j) ||
                    (i >= size1_))
                    break;
                ++ i;
                it = data ().lower_bound (functor_type::element (i, size1_,
j, size2_));
            }
            if (rank == 1 && index2 != j) {
                i = size1_; <---------- Add this line
                rank = 0;
            }
            return const_iterator1 (*this, rank, i, j, it);
        }
        iterator1 find1 (int rank, size_type i, size_type j) {
            iterator_type it (data ().lower_bound (functor_type::element (i,
size1_, j, size2_)));
            iterator_type it_end (data ().end ());
            size_type index2 = size_type (-1);
            while (it != it_end) {
                index2 = functor_type::index2 ((*it).first, size1_, size2_);
                if ((rank == 0 && index2 >= j) ||
                    (rank == 1 && index2 == j) ||
                    (i >= size1_))
                    break;
                ++ i;
                it = data ().lower_bound (functor_type::element (i, size1_,
j, size2_));
            }
            if (rank == 1 && index2 != j) {
                i = size1_; <---------- Add this line
                rank = 0;
            }
            return iterator1 (*this, rank, i, j, it);
        }
        const_iterator2 find2 (int rank, size_type i, size_type j) const {
            const_iterator_type it (data ().lower_bound
(functor_type::element (i, size1_, j, size2_)));
            const_iterator_type it_end (data ().end ());
            size_type index1 = size_type (-1);
            while (it != it_end) {
                index1 = functor_type::index1 ((*it).first, size1_, size2_);
                if ((rank == 0 && index1 >= i) ||
                    (rank == 1 && index1 == i) ||
                    (j >= size2_)) {
                    break;
                }
                ++ j;
                it = data ().lower_bound (functor_type::element (i, size1_,
j, size2_));
            }
            if (rank == 1 && index1 != i) {
                j = size2_; <---------- Add this line
                rank = 0;
            }
            return const_iterator2 (*this, rank, i, j, it);
        }
        iterator2 find2 (int rank, size_type i, size_type j) {
            iterator_type it (data ().lower_bound (functor_type::element (i,
size1_, j, size2_)));
            iterator_type it_end (data ().end ());
            size_type index1 = size_type (-1);
            while (it != it_end) {
                index1 = functor_type::index1 ((*it).first, size1_, size2_);
                if ((rank == 0 && index1 >= i) ||
                    (rank == 1 && index1 == i) ||
                    (j >= size2_)) {
                    break;
                }
                ++ j;
                it = data ().lower_bound (functor_type::element (i, size1_,
j, size2_));
            }
            if (rank == 1 && index1 != i) {
                j = size2_; <---------- Add this line
                rank = 0;
            }
            return iterator2 (*this, rank, i, j, it);
        }

Thanks for your feedback

Joerg

P.S.: I've run my basic tests against this change.

P.P.S.: I'm not sure, whether I should check in our current version to boost
CVS during review.