From 02c8b6af821fddf1976141448977dfb12c71ec30 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Wed, 27 Oct 2010 15:13:03 +0200
Subject: [PATCH] fix sparse rankUpdate and triangularView iterator

---
 Eigen/src/Sparse/SparseSelfAdjointView.h | 10 +++++-----
 Eigen/src/Sparse/SparseTriangularView.h  | 14 ++++++++------
 unsupported/test/sparse_llt.cpp          | 13 +++++++++++--
 3 files changed, 24 insertions(+), 13 deletions(-)

diff --git a/Eigen/src/Sparse/SparseSelfAdjointView.h b/Eigen/src/Sparse/SparseSelfAdjointView.h
index 769a05cf7..2a7101082 100644
--- a/Eigen/src/Sparse/SparseSelfAdjointView.h
+++ b/Eigen/src/Sparse/SparseSelfAdjointView.h
@@ -92,7 +92,7 @@ template<typename MatrixType, unsigned int UpLo> class SparseSelfAdjointView
       * call this function with u.adjoint().
       */
     template<typename DerivedU>
-    SparseSelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, Scalar alpha = Scalar(1));
+    SparseSelfAdjointView& rankUpdate(const SparseMatrixBase<DerivedU>& u, Scalar alpha = Scalar(1));
 
     // const SparseLLT<PlainObject, UpLo> llt() const;
     // const SparseLDLT<PlainObject, UpLo> ldlt() const;
@@ -127,15 +127,15 @@ SparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView(
 template<typename MatrixType, unsigned int UpLo>
 template<typename DerivedU>
 SparseSelfAdjointView<MatrixType,UpLo>&
-SparseSelfAdjointView<MatrixType,UpLo>::rankUpdate(const MatrixBase<DerivedU>& u, Scalar alpha)
+SparseSelfAdjointView<MatrixType,UpLo>::rankUpdate(const SparseMatrixBase<DerivedU>& u, Scalar alpha)
 {
   SparseMatrix<Scalar,MatrixType::Flags&RowMajorBit?RowMajor:ColMajor> tmp = u * u.adjoint();
   if(alpha==Scalar(0))
-    m_matrix = tmp.template triangularView<UpLo>();
+    m_matrix.const_cast_derived() = tmp.template triangularView<UpLo>();
   else
-    m_matrix += alpha * tmp.template triangularView<UpLo>();
+    m_matrix.const_cast_derived() /*+*/= alpha * tmp.template triangularView<UpLo>();
 
-  return this;
+  return *this;
 }
 
 /***************************************************************************
diff --git a/Eigen/src/Sparse/SparseTriangularView.h b/Eigen/src/Sparse/SparseTriangularView.h
index 8d0256599..319eaf066 100644
--- a/Eigen/src/Sparse/SparseTriangularView.h
+++ b/Eigen/src/Sparse/SparseTriangularView.h
@@ -26,11 +26,13 @@
 #define EIGEN_SPARSE_TRIANGULARVIEW_H
 
 namespace internal {
+  
 template<typename MatrixType, int Mode>
 struct traits<SparseTriangularView<MatrixType,Mode> >
 : public traits<MatrixType>
 {};
-}
+
+} // namespace internal
 
 template<typename MatrixType, int Mode> class SparseTriangularView
   : public SparseMatrixBase<SparseTriangularView<MatrixType,Mode> >
@@ -38,13 +40,13 @@ template<typename MatrixType, int Mode> class SparseTriangularView
     enum { SkipFirst = (Mode==Lower && !(MatrixType::Flags&RowMajorBit))
                     || (Mode==Upper &&  (MatrixType::Flags&RowMajorBit)) };
   public:
+    
+    EIGEN_SPARSE_PUBLIC_INTERFACE(SparseTriangularView)
 
     class InnerIterator;
-    typedef typename MatrixType::Scalar Scalar;
-    typedef typename MatrixType::Index Index;
 
-    inline Index rows() { return m_matrix.rows(); }
-    inline Index cols() { return m_matrix.cols(); }
+    inline Index rows() const { return m_matrix.rows(); }
+    inline Index cols() const { return m_matrix.cols(); }
 
     typedef typename internal::conditional<internal::must_nest_by_value<MatrixType>::ret,
         MatrixType, const MatrixType&>::type MatrixTypeNested;
@@ -83,7 +85,7 @@ class SparseTriangularView<MatrixType,Mode>::InnerIterator : public MatrixType::
 
     EIGEN_STRONG_INLINE operator bool() const
     {
-      return SkipFirst ? Base::operator bool() : (Base::operator bool() && this->index() < this->outer());
+      return SkipFirst ? Base::operator bool() : (Base::operator bool() && this->index() <= this->outer());
     }
 };
 
diff --git a/unsupported/test/sparse_llt.cpp b/unsupported/test/sparse_llt.cpp
index 7b1cad4c5..2fcd8b279 100644
--- a/unsupported/test/sparse_llt.cpp
+++ b/unsupported/test/sparse_llt.cpp
@@ -77,8 +77,11 @@ template<typename Scalar> void sparse_llt(int rows, int cols)
     // new API
     {
       // Cholmod, as configured in CholmodSupport.h, only supports self-adjoint matrices
-      SparseMatrix<Scalar> m3 = m2.adjoint()*m2;
-      DenseMatrix refMat3 = refMat2.adjoint()*refMat2;
+      SparseMatrix<Scalar> m3 = m2 * m2.adjoint(), m3_lo(rows,rows), m3_up(rows,rows);
+      DenseMatrix refMat3 = refMat2 * refMat2.adjoint();
+      
+      m3_lo.template selfadjointView<Lower>().rankUpdate(m2,0);
+      m3_up.template selfadjointView<Upper>().rankUpdate(m2,0);
       
       // with a single vector as the rhs
       ref_x = refMat3.template selfadjointView<Lower>().llt().solve(b);
@@ -89,6 +92,12 @@ template<typename Scalar> void sparse_llt(int rows, int cols)
       x = CholmodDecomposition<SparseMatrix<Scalar>, Upper>(m3).solve(b);
       VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod solve, single dense rhs");
       
+      x = CholmodDecomposition<SparseMatrix<Scalar>, Lower>(m3_lo).solve(b);
+      VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod solve, single dense rhs");
+      
+      x = CholmodDecomposition<SparseMatrix<Scalar>, Upper>(m3_up).solve(b);
+      VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod solve, single dense rhs");
+      
       
       // with multiple rhs
       ref_X = refMat3.template selfadjointView<Lower>().llt().solve(B);