[eigen] [Patch] Bug in Sparse module

[David Luitz <tux008@xxxxxxxxxxxxxx>]

Hi all,

There is a bug in the Sparse Module, in 
Eigen/src/Sparse/SparseDenseProduct.h class SparseTimeDenseProduct.

The problem is, that j runs from 0 to lhs.outerSize(), and coeff(j,0) 
may be out of range if lhs is Row major. This results in a segmentation 
fault, which is illustrated by the piece of code attached in the file 
sparse2.cc. (Change the commented line to compare both cases.)

I have written a patch that solves this problem by changing the order of 
the existing code. The new code also eliminates the temporary variables 
dest_j and rhs_j in the cases where they are not needed.

The patch should work against mercurial revision 3604:537b5d33b440

Greetings,
David J. Luitz

#include<iostream>
#include <Eigen/Eigen>
int main(){
	Eigen::SparseMatrix<double, Eigen::RowMajor> m1(10000, 2);// segfaults
//	Eigen::SparseMatrix<double, Eigen::ColMajor> m1(10000, 2);// works
	Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> m2(2, 2);
	m1.setZero();
	m2.setZero(2,2);

	Eigen::Matrix<double,Eigen::Dynamic, Eigen::Dynamic> res( m1*m2 );

	std::cout << res(0,0) << std::endl;
}

diff -r 537b5d33b440 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:41:17 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;