[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
- To: eigen@xxxxxxxxxxxxxxxxxxx
- Subject: [eigen] [Sparse] bugfix
- From: David Luitz <tux008@xxxxxxxxxxxxxx>
- Date: Thu, 30 Dec 2010 15:25:42 +0100
- Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=googlemail.com; s=gamma; h=domainkey-signature:received:received:message-id:date:from :user-agent:mime-version:to:subject:content-type; bh=UwAtmeWZ7jkoOr8ioeUyncU/VQ8zfg4YYDL6LhEhkAw=; b=b2t5Hio0SUhWUdCuSXM97mUHumqdysD2ffp1ejzCBDPE+xiSrILWmEfkbXckNh2NES Sbwb8iUvm8Qp6m9hstwXiWWmXtS4rSU9HDKJubDFZfC7wALsSI3GWowsmn30xn9EH+pL RgR8zdv14FGmHMGSii34fkhy6bct/igndliDM=
- Domainkey-signature: a=rsa-sha1; c=nofws; d=googlemail.com; s=gamma; h=message-id:date:from:user-agent:mime-version:to:subject :content-type; b=XheAP1pLGAqzDO/lbqi3XuvZiveVF8a8W1Z17haa8J3f3RE3vKfinpVUGr1F2SCCTz HI0pQ+MTwDATbEJFepsLogJ4hVuXInf4567FEC8sBa5epmjW/8vDxGUvzcaA7kedh/II btmuIxQ2Gh3LjUL8Iaa71ckyYlKmOJEZhqnFQ=
Me again,
I have another fix for a bug in the sparse module, same place as before
in the file Eigen/src/Sparse/SparseDenseProduct.h.
In the class SparseTimeDenseProduct, we use a Block constructor, which
does not compile in certain cases.
Block<Dest,1,Dest::ColsAtCompileTime> dest_j(dest.row(LhsIsRowMajor ? j
: 0))
I changed it to
Block<Dest,1,Dest::ColsAtCompileTime> dest_j(dest, LhsIsRowMajor ? j : 0);
and this seems to work now. I added a unit test sparse_product_4, which
fails to compile with the old code and works fine with the new code. The
test also checks for the Sparse bug I reported before and produces a
segfault (at least on my machine) for the old code, while it works with
my patch.
I commited all these changes including my previous patch to mercurial
and exported them using hg export.
You will find 3 commits in the appended patch file.
Greetings,
David Luitz
exporting patches:
# HG changeset patch
# User David J. Luitz <tux008@xxxxxxxxxxxxxx>
# Date 1293718583 -3600
# Node ID c6ef89f7282eeb83035fa946747d293ffeb3513a
# Parent 537b5d33b440cc42188476a1e7d1894a45cd4910
# Parent 8937cfef58d27ce57e307c5531d4307441f91107
[Sparse] Added regression tests for the two bugfixes, the code passes all sparse_product tests
diff -r 537b5d33b440 -r c6ef89f7282e Eigen/src/Sparse/SparseDenseProduct.h
--- a/Eigen/src/Sparse/SparseDenseProduct.h Tue Dec 28 13:46:39 2010 +0100
+++ b/Eigen/src/Sparse/SparseDenseProduct.h Thu Dec 30 15:16:23 2010 +0100
@@ -169,14 +169,31 @@
enum { LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit };
for(Index j=0; j<m_lhs.outerSize(); ++j)
{
- typename Rhs::Scalar rhs_j = alpha * m_rhs.coeff(j,0);
- Block<Dest,1,Dest::ColsAtCompileTime> dest_j(dest.row(LhsIsRowMajor ? j : 0));
- for(LhsInnerIterator it(m_lhs,j); it ;++it)
- {
- if(LhsIsRowMajor) dest_j += (alpha*it.value()) * m_rhs.row(it.index());
- else if(Rhs::ColsAtCompileTime==1) dest.coeffRef(it.index()) += it.value() * rhs_j;
- else dest.row(it.index()) += (alpha*it.value()) * m_rhs.row(j);
- }
+ if(LhsIsRowMajor)
+ {
+// Block<Dest,1,Dest::ColsAtCompileTime> dest_j(dest.row(LhsIsRowMajor ? j : 0)); // this does not work in all cases. Why?
+ Block<Dest,1,Dest::ColsAtCompileTime> dest_j(dest, LhsIsRowMajor ? j : 0);
+ for(LhsInnerIterator it(m_lhs,j); it ;++it)
+ {
+ dest_j += (alpha*it.value()) * m_rhs.row(it.index());
+ }
+ }
+ else if(Rhs::ColsAtCompileTime==1)
+ {
+ typename Rhs::Scalar rhs_j = alpha * m_rhs.coeff(j,0);
+ for(LhsInnerIterator it(m_lhs,j); it ;++it)
+ {
+ dest.coeffRef(it.index()) += it.value() * rhs_j;
+ }
+ }
+ else
+ {
+ for(LhsInnerIterator it(m_lhs,j); it ;++it)
+ {
+ dest.row(it.index()) += (alpha*it.value()) * m_rhs.row(j);
+ }
+ }
+
}
}
@@ -187,8 +204,8 @@
// dense = dense * sparse
namespace internal {
-template<typename Lhs, typename Rhs>
-struct traits<DenseTimeSparseProduct<Lhs,Rhs> >
+ template<typename Lhs, typename Rhs>
+ struct traits<DenseTimeSparseProduct<Lhs,Rhs> >
: traits<ProductBase<DenseTimeSparseProduct<Lhs,Rhs>, Lhs, Rhs> >
{
typedef Dense StorageKind;
diff -r 537b5d33b440 -r c6ef89f7282e test/sparse_product.cpp
--- a/test/sparse_product.cpp Tue Dec 28 13:46:39 2010 +0100
+++ b/test/sparse_product.cpp Thu Dec 30 15:16:23 2010 +0100
@@ -137,6 +137,31 @@
}
}
+// New test for Bug in SparseTimeDenseProduct
+template<typename SparseMatrixType, typename DenseMatrixType> void sparse_product_regression_test()
+{
+ // This code does not compile with afflicted versions of the bug
+/* SparseMatrixType sm1(3,2);
+ DenseMatrixType m2(2,2);
+ sm1.setZero();
+ m2.setZero();
+
+ DenseMatrixType m3 = sm1*m2;
+ */
+
+
+ // This code produces a segfault with afflicted versions of another SparseTimeDenseProduct
+ // bug
+
+ SparseMatrixType sm2(20000,2);
+ DenseMatrixType m3(2,2);
+ sm2.setZero();
+ m3.setZero();
+ DenseMatrixType m4(sm2*m3);
+
+ VERIFY_IS_APPROX( m4(0,0), 0.0 );
+}
+
void test_sparse_product()
{
for(int i = 0; i < g_repeat; i++) {
@@ -145,5 +170,7 @@
CALL_SUBTEST_1( sparse_product(SparseMatrix<double>(33, 33)) );
CALL_SUBTEST_3( sparse_product(DynamicSparseMatrix<double>(8, 8)) );
+
+ CALL_SUBTEST_4( (sparse_product_regression_test<SparseMatrix<double,RowMajor>, Matrix<double, Dynamic, Dynamic, RowMajor> >()) );
}
}
# HG changeset patch
# User David J. Luitz <tux008@xxxxxxxxxxxxxx>
# Date 1293702327 -3600
# Branch experimental
# Node ID 8937cfef58d27ce57e307c5531d4307441f91107
# Parent 8218e6577d32463d3902e4d66431a5833ef12875
[Sparse] Compile fix in SparseDenseProduct.
+ Apparently, Block<> b( dest.row(j) ) does not work in all cases. I changed this line to use
another constructor.
diff -r 8218e6577d32 -r 8937cfef58d2 Eigen/src/Sparse/SparseDenseProduct.h
--- a/Eigen/src/Sparse/SparseDenseProduct.h Wed Dec 29 00:35:00 2010 +0100
+++ b/Eigen/src/Sparse/SparseDenseProduct.h Thu Dec 30 10:45:27 2010 +0100
@@ -171,7 +171,8 @@
{
if(LhsIsRowMajor)
{
- Block<Dest,1,Dest::ColsAtCompileTime> dest_j(dest.row(LhsIsRowMajor ? j : 0));
+// Block<Dest,1,Dest::ColsAtCompileTime> dest_j(dest.row(LhsIsRowMajor ? j : 0)); // this does not work in all cases. Why?
+ Block<Dest,1,Dest::ColsAtCompileTime> dest_j(dest, LhsIsRowMajor ? j : 0);
for(LhsInnerIterator it(m_lhs,j); it ;++it)
{
dest_j += (alpha*it.value()) * m_rhs.row(it.index());
# HG changeset patch
# User David J. Luitz <tux008@xxxxxxxxxxxxxx>
# Date 1293579300 -3600
# Branch experimental
# Node ID 8218e6577d32463d3902e4d66431a5833ef12875
# Parent 537b5d33b440cc42188476a1e7d1894a45cd4910
[Sparse] bugfix for SparseTimeDenseProduct
diff -r 537b5d33b440 -r 8218e6577d32 Eigen/src/Sparse/SparseDenseProduct.h
--- a/Eigen/src/Sparse/SparseDenseProduct.h Tue Dec 28 13:46:39 2010 +0100
+++ b/Eigen/src/Sparse/SparseDenseProduct.h Wed Dec 29 00:35:00 2010 +0100
@@ -169,14 +169,30 @@
enum { LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit };
for(Index j=0; j<m_lhs.outerSize(); ++j)
{
- typename Rhs::Scalar rhs_j = alpha * m_rhs.coeff(j,0);
- Block<Dest,1,Dest::ColsAtCompileTime> dest_j(dest.row(LhsIsRowMajor ? j : 0));
- for(LhsInnerIterator it(m_lhs,j); it ;++it)
- {
- if(LhsIsRowMajor) dest_j += (alpha*it.value()) * m_rhs.row(it.index());
- else if(Rhs::ColsAtCompileTime==1) dest.coeffRef(it.index()) += it.value() * rhs_j;
- else dest.row(it.index()) += (alpha*it.value()) * m_rhs.row(j);
- }
+ if(LhsIsRowMajor)
+ {
+ Block<Dest,1,Dest::ColsAtCompileTime> dest_j(dest.row(LhsIsRowMajor ? j : 0));
+ for(LhsInnerIterator it(m_lhs,j); it ;++it)
+ {
+ dest_j += (alpha*it.value()) * m_rhs.row(it.index());
+ }
+ }
+ else if(Rhs::ColsAtCompileTime==1)
+ {
+ typename Rhs::Scalar rhs_j = alpha * m_rhs.coeff(j,0);
+ for(LhsInnerIterator it(m_lhs,j); it ;++it)
+ {
+ dest.coeffRef(it.index()) += it.value() * rhs_j;
+ }
+ }
+ else
+ {
+ for(LhsInnerIterator it(m_lhs,j); it ;++it)
+ {
+ dest.row(it.index()) += (alpha*it.value()) * m_rhs.row(j);
+ }
+ }
+
}
}
@@ -187,8 +203,8 @@
// dense = dense * sparse
namespace internal {
-template<typename Lhs, typename Rhs>
-struct traits<DenseTimeSparseProduct<Lhs,Rhs> >
+ template<typename Lhs, typename Rhs>
+ struct traits<DenseTimeSparseProduct<Lhs,Rhs> >
: traits<ProductBase<DenseTimeSparseProduct<Lhs,Rhs>, Lhs, Rhs> >
{
typedef Dense StorageKind;