From 0a537cb2d87ada8206ec2271fb9f2904a18ccfce Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Fri, 12 Feb 2016 15:58:31 +0100
Subject: [PATCH] bug #901: fix triangular-view with unit diagonal of sparse
 rectangular matrices.

---
 Eigen/src/SparseCore/SparseTriangularView.h | 21 ++++++++++++---------
 test/sparse_basic.cpp                       |  5 ++---
 2 files changed, 14 insertions(+), 12 deletions(-)

diff --git a/Eigen/src/SparseCore/SparseTriangularView.h b/Eigen/src/SparseCore/SparseTriangularView.h
index 7c718e4e1..2c6aedaf9 100644
--- a/Eigen/src/SparseCore/SparseTriangularView.h
+++ b/Eigen/src/SparseCore/SparseTriangularView.h
@@ -70,20 +70,20 @@ class TriangularViewImpl<MatrixType,Mode,Sparse>::InnerIterator : public MatrixT
   public:
 
     EIGEN_STRONG_INLINE InnerIterator(const TriangularViewImpl& view, Index outer)
-      : Base(view.derived().nestedExpression(), outer), m_returnOne(false)
+      : Base(view.derived().nestedExpression(), outer), m_returnOne(false), m_containsDiag(Base::outer()<view.innerSize())
     {
       if(SkipFirst)
       {
         while((*this) && ((HasUnitDiag||SkipDiag)  ? this->index()<=outer : this->index()<outer))
           Base::operator++();
         if(HasUnitDiag)
-          m_returnOne = true;
+          m_returnOne = m_containsDiag;
       }
       else if(HasUnitDiag && ((!Base::operator bool()) || Base::index()>=Base::outer()))
       {
         if((!SkipFirst) && Base::operator bool())
           Base::operator++();
-        m_returnOne = true;
+        m_returnOne = m_containsDiag;
       }
     }
 
@@ -98,7 +98,7 @@ class TriangularViewImpl<MatrixType,Mode,Sparse>::InnerIterator : public MatrixT
         {
           if((!SkipFirst) && Base::operator bool())
             Base::operator++();
-          m_returnOne = true;
+          m_returnOne = m_containsDiag;
         }
       }
       return *this;
@@ -130,6 +130,7 @@ class TriangularViewImpl<MatrixType,Mode,Sparse>::InnerIterator : public MatrixT
     }
   protected:
     bool m_returnOne;
+    bool m_containsDiag;
 };
 
 template<typename MatrixType, unsigned int Mode>
@@ -193,7 +194,7 @@ public:
     Flags = XprType::Flags
   };
     
-  explicit unary_evaluator(const XprType &xpr) : m_argImpl(xpr.nestedExpression()) {}
+  explicit unary_evaluator(const XprType &xpr) : m_argImpl(xpr.nestedExpression()), m_arg(xpr.nestedExpression()) {}
   
   inline Index nonZerosEstimate() const {
     return m_argImpl.nonZerosEstimate();
@@ -205,20 +206,20 @@ public:
     public:
 
       EIGEN_STRONG_INLINE InnerIterator(const unary_evaluator& xprEval, Index outer)
-        : Base(xprEval.m_argImpl,outer), m_returnOne(false)
+        : Base(xprEval.m_argImpl,outer), m_returnOne(false), m_containsDiag(Base::outer()<xprEval.m_arg.innerSize())
       {
         if(SkipFirst)
         {
           while((*this) && ((HasUnitDiag||SkipDiag)  ? this->index()<=outer : this->index()<outer))
             Base::operator++();
           if(HasUnitDiag)
-            m_returnOne = true;
+            m_returnOne = m_containsDiag;
         }
         else if(HasUnitDiag && ((!Base::operator bool()) || Base::index()>=Base::outer()))
         {
           if((!SkipFirst) && Base::operator bool())
             Base::operator++();
-          m_returnOne = true; // FIXME check innerSize()>outer();
+          m_returnOne = m_containsDiag;
         }
       }
 
@@ -233,7 +234,7 @@ public:
           {
             if((!SkipFirst) && Base::operator bool())
               Base::operator++();
-            m_returnOne = true; // FIXME check innerSize()>outer();
+            m_returnOne = m_containsDiag;
           }
         }
         return *this;
@@ -266,12 +267,14 @@ public:
 
     protected:
       bool m_returnOne;
+      bool m_containsDiag;
     private:
       Scalar& valueRef();
   };
   
 protected:
   evaluator<ArgType> m_argImpl;
+  const ArgType& m_arg;
 };
 
 } // end namespace internal
diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp
index 0a06c828b..cb8ebaedf 100644
--- a/test/sparse_basic.cpp
+++ b/test/sparse_basic.cpp
@@ -21,8 +21,8 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re
   
   const Index rows = ref.rows();
   const Index cols = ref.cols();
-  const Index inner = ref.innerSize();
-  const Index outer = ref.outerSize();
+  //const Index inner = ref.innerSize();
+  //const Index outer = ref.outerSize();
 
   typedef typename SparseMatrixType::Scalar Scalar;
   enum { Flags = SparseMatrixType::Flags };
@@ -327,7 +327,6 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re
     m3 = m2.template triangularView<Upper>();
     VERIFY_IS_APPROX(m3, refMat3);
 
-    if(inner>=outer) // FIXME this should be implemented for outer>inner as well
     {
       refMat3 = refMat2.template triangularView<UnitUpper>();
       m3 = m2.template triangularView<UnitUpper>();