# Boost :

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);
}