From 296d24be4dd3c700089d1d9182a843c60d68019c Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Wed, 25 Jan 2017 17:39:01 +0100
Subject: [PATCH] bug #1381: fix sparse.diagonal() used as a rvalue. The
 problem was that is "sparse" is not const, then sparse.diagonal() must have
 the LValueBit flag meaning that sparse.diagonal().coeff(i) must returns a
 const reference, const Scalar&. However, sparse::coeff() cannot returns a
 reference for a non-existing zero coefficient. The trick is to return a
 reference to a local member of evaluator<SparseMatrix>.

---
 Eigen/src/Core/CoreEvaluators.h             |  4 +--
 Eigen/src/Core/util/XprHelper.h             |  2 +-
 Eigen/src/SparseCore/SparseCompressedBase.h | 38 +++++++++++++++------
 test/sparse_basic.cpp                       |  4 +++
 4 files changed, 33 insertions(+), 15 deletions(-)

diff --git a/Eigen/src/Core/CoreEvaluators.h b/Eigen/src/Core/CoreEvaluators.h
index 24fc7835b..412f5a661 100644
--- a/Eigen/src/Core/CoreEvaluators.h
+++ b/Eigen/src/Core/CoreEvaluators.h
@@ -1613,9 +1613,7 @@ struct evaluator<Diagonal<ArgType, DiagIndex> >
   { }
  
   typedef typename XprType::Scalar Scalar;
-  // FIXME having to check whether ArgType is sparse here i not very nice.
-  typedef typename internal::conditional<!internal::is_same<typename ArgType::StorageKind,Sparse>::value,
-                                         typename XprType::CoeffReturnType,Scalar>::type CoeffReturnType;
+  typedef typename XprType::CoeffReturnType CoeffReturnType;
 
   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
   CoeffReturnType coeff(Index row, Index) const
diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h
index b94dd61ff..d5e565cb8 100644
--- a/Eigen/src/Core/util/XprHelper.h
+++ b/Eigen/src/Core/util/XprHelper.h
@@ -638,7 +638,7 @@ struct plain_constant_type
 template<typename ExpressionType>
 struct is_lvalue
 {
-  enum { value = !bool(is_const<ExpressionType>::value) &&
+  enum { value = (!bool(is_const<ExpressionType>::value)) &&
                  bool(traits<ExpressionType>::Flags & LvalueBit) };
 };
 
diff --git a/Eigen/src/SparseCore/SparseCompressedBase.h b/Eigen/src/SparseCore/SparseCompressedBase.h
index d32a8df40..e0b3c22b6 100644
--- a/Eigen/src/SparseCore/SparseCompressedBase.h
+++ b/Eigen/src/SparseCore/SparseCompressedBase.h
@@ -295,11 +295,11 @@ struct evaluator<SparseCompressedBase<Derived> >
     Flags = Derived::Flags
   };
   
-  evaluator() : m_matrix(0)
+  evaluator() : m_matrix(0), m_zero(0)
   {
     EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost);
   }
-  explicit evaluator(const Derived &mat) : m_matrix(&mat)
+  explicit evaluator(const Derived &mat) : m_matrix(&mat), m_zero(0)
   {
     EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost);
   }
@@ -312,26 +312,42 @@ struct evaluator<SparseCompressedBase<Derived> >
   operator const Derived&() const { return *m_matrix; }
   
   typedef typename DenseCoeffsBase<Derived,ReadOnlyAccessors>::CoeffReturnType CoeffReturnType;
-  Scalar coeff(Index row, Index col) const
-  { return m_matrix->coeff(row,col); }
-  
+  const Scalar& coeff(Index row, Index col) const
+  {
+    Index p = find(row,col);
+
+    if(p==Dynamic)
+      return m_zero;
+    else
+      return m_matrix->const_cast_derived().valuePtr()[p];
+  }
+
   Scalar& coeffRef(Index row, Index col)
+  {
+    Index p = find(row,col);
+    eigen_assert(p!=Dynamic && "written coefficient does not exist");
+    return m_matrix->const_cast_derived().valuePtr()[p];
+  }
+
+protected:
+
+  Index find(Index row, Index col) const
   {
     eigen_internal_assert(row>=0 && row<m_matrix->rows() && col>=0 && col<m_matrix->cols());
-      
+
     const Index outer = Derived::IsRowMajor ? row : col;
     const Index inner = Derived::IsRowMajor ? col : row;
 
     Index start = m_matrix->outerIndexPtr()[outer];
     Index end = m_matrix->isCompressed() ? m_matrix->outerIndexPtr()[outer+1] : m_matrix->outerIndexPtr()[outer] + m_matrix->innerNonZeroPtr()[outer];
-    eigen_assert(end>start && "you are using a non finalized sparse matrix or written coefficient does not exist");
-    const Index p =   std::lower_bound(m_matrix->innerIndexPtr()+start, m_matrix->innerIndexPtr()+end,inner)
-                    - m_matrix->innerIndexPtr();
-    eigen_assert((p<end) && (m_matrix->innerIndexPtr()[p]==inner) && "written coefficient does not exist");
-    return m_matrix->const_cast_derived().valuePtr()[p];
+    eigen_assert(end>=start && "you are using a non finalized sparse matrix or written coefficient does not exist");
+    const Index p = std::lower_bound(m_matrix->innerIndexPtr()+start, m_matrix->innerIndexPtr()+end,inner) - m_matrix->innerIndexPtr();
+
+    return ((p<end) && (m_matrix->innerIndexPtr()[p]==inner)) ? p : Dynamic;
   }
 
   const Derived *m_matrix;
+  const Scalar m_zero;
 };
 
 }
diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp
index 208a66f12..91b7cb335 100644
--- a/test/sparse_basic.cpp
+++ b/test/sparse_basic.cpp
@@ -485,6 +485,10 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re
     SparseMatrixType m2(rows, cols);
     initSparse<Scalar>(density, refMat2, m2);
     VERIFY_IS_APPROX(m2.diagonal(), refMat2.diagonal().eval());
+    DenseVector d = m2.diagonal();
+    VERIFY_IS_APPROX(d, refMat2.diagonal().eval());
+    d = m2.diagonal().array();
+    VERIFY_IS_APPROX(d, refMat2.diagonal().eval());
     VERIFY_IS_APPROX(const_cast<const SparseMatrixType&>(m2).diagonal(), refMat2.diagonal().eval());
     
     initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag);