new feature: copy from a sparse selfadjoint view to a full sparse matrix

5 files changed, 104 insertions, 4 deletions

diff --git a/Eigen/src/Sparse/DynamicSparseMatrix.h b/Eigen/src/Sparse/DynamicSparseMatrix.h
index 88ad85f84..49a78fd26 100644
--- a/Eigen/src/Sparse/DynamicSparseMatrix.h
+++ b/Eigen/src/Sparse/DynamicSparseMatrix.h
@@ -108,7 +108,7 @@ class DynamicSparseMatrix
 
     /** \returns a reference to the coefficient value at given position \a row, \a col
       * This operation involes a log(rho*outer_size) binary search. If the coefficient does not
-      * exist yet, then a sorted insertion Indexo a sequential buffer is performed.
+      * exist yet, then a sorted insertion into a sequential buffer is performed.
       */
     inline Scalar& coeffRef(Index row, Index col)
     {
diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h
index 38e810725..4fa5d5e33 100644
--- a/Eigen/src/Sparse/SparseMatrix.h
+++ b/Eigen/src/Sparse/SparseMatrix.h
@@ -69,6 +69,7 @@ class SparseMatrix
 {
   public:
     EIGEN_SPARSE_PUBLIC_INTERFACE(SparseMatrix)
+    using Base::operator=;
     EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(SparseMatrix, +=)
     EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(SparseMatrix, -=)
     // FIXME: why are these operator already alvailable ???
diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h
index ba6f64ab7..bc16be955 100644
--- a/Eigen/src/Sparse/SparseMatrixBase.h
+++ b/Eigen/src/Sparse/SparseMatrixBase.h
@@ -46,6 +46,16 @@ template<typename Derived> class SparseMatrixBase : public EigenBase<Derived>
     typedef typename internal::traits<Derived>::Index Index;
 
     typedef SparseMatrixBase StorageBaseType;
+    typedef EigenBase<Derived> Base;
+    
+    template<typename OtherDerived>
+    Derived& operator=(const EigenBase<OtherDerived> &other)
+    {
+      other.derived().evalTo(derived());
+      return derived();
+    }
+    
+//     using Base::operator=;
 
     enum {
 
@@ -173,7 +183,7 @@ template<typename Derived> class SparseMatrixBase : public EigenBase<Derived>
     Derived& markAsRValue() { m_isRValue = true; return derived(); }
 
     SparseMatrixBase() : m_isRValue(false) { /* TODO check flags */ }
-
+    
     inline Derived& operator=(const Derived& other)
     {
 //       std::cout << "Derived& operator=(const Derived& other)\n";
diff --git a/Eigen/src/Sparse/SparseSelfAdjointView.h b/Eigen/src/Sparse/SparseSelfAdjointView.h
index 2a7101082..b4ece8dd4 100644
--- a/Eigen/src/Sparse/SparseSelfAdjointView.h
+++ b/Eigen/src/Sparse/SparseSelfAdjointView.h
@@ -45,12 +45,24 @@ class SparseSelfAdjointTimeDenseProduct;
 template<typename Lhs, typename Rhs, int UpLo>
 class DenseTimeSparseSelfAdjointProduct;
 
+namespace internal {
+  
+template<typename MatrixType, unsigned int UpLo>
+struct traits<SparseSelfAdjointView<MatrixType,UpLo> > : traits<MatrixType> {
+};
+
+}
+
 template<typename MatrixType, unsigned int UpLo> class SparseSelfAdjointView
+  : public EigenBase<SparseSelfAdjointView<MatrixType,UpLo> >
 {
   public:
 
     typedef typename MatrixType::Scalar Scalar;
     typedef typename MatrixType::Index Index;
+    typedef Matrix<Index,Dynamic,1> VectorI;
+    typedef typename MatrixType::Nested MatrixTypeNested;
+    typedef typename internal::remove_all<MatrixTypeNested>::type _MatrixTypeNested;
 
     inline SparseSelfAdjointView(const MatrixType& matrix) : m_matrix(matrix)
     {
@@ -77,7 +89,7 @@ template<typename MatrixType, unsigned int UpLo> class SparseSelfAdjointView
     DenseTimeSparseSelfAdjointProduct<OtherDerived,MatrixType,UpLo>
     operator*(const MatrixBase<OtherDerived>& lhs, const SparseSelfAdjointView& rhs)
     {
-      return DenseTimeSparseSelfAdjointProduct<OtherDerived,MatrixType,UpLo>(lhs.derived(), rhs.m_matrix);
+      return DenseTimeSparseSelfAdjointProduct<OtherDerived,_MatrixTypeNested,UpLo>(lhs.derived(), rhs.m_matrix);
     }
 
     /** Perform a symmetric rank K update of the selfadjoint matrix \c *this:
@@ -93,6 +105,70 @@ template<typename MatrixType, unsigned int UpLo> class SparseSelfAdjointView
       */
     template<typename DerivedU>
     SparseSelfAdjointView& rankUpdate(const SparseMatrixBase<DerivedU>& u, Scalar alpha = Scalar(1));
+    
+    /** \internal triggered by sparse_matrix = SparseSelfadjointView; */
+    template<typename DestScalar> void evalTo(SparseMatrix<DestScalar>& _dest) const
+    {
+      typedef SparseMatrix<DestScalar> Dest;
+      Dest& dest(_dest.derived());
+      enum {
+        StorageOrderMatch = Dest::IsRowMajor == _MatrixTypeNested::IsRowMajor
+      };
+      VectorI count = Dest::IsRowMajor ? m_countPerRow : m_countPerCol;
+      Index size = m_matrix.rows();
+      count.resize(size);
+      count.setZero();
+      dest.resize(size,size);
+      for(Index j = 0; j<size; ++j)
+      {
+        for(typename _MatrixTypeNested::InnerIterator it(m_matrix,j); it; ++it)
+        {
+          Index i = it.index();
+          if(i==j)
+            count[i]++;
+          else if((UpLo==Lower && i>j) || (UpLo==Upper && i<j))
+          {
+            count[i]++;
+            count[j]++;
+          }
+        }
+      }
+      Index nnz = count.sum();
+      
+      // reserve space
+      dest.reserve(nnz);
+      dest._outerIndexPtr()[0] = 0;
+      for(Index j=0; j<size; ++j)
+        dest._outerIndexPtr()[j+1] = dest._outerIndexPtr()[j] + count[j];
+      for(Index j=0; j<size; ++j)
+        count[j] = dest._outerIndexPtr()[j];
+      
+      // copy data
+      for(Index j = 0; j<size; ++j)
+      {
+        for(typename _MatrixTypeNested::InnerIterator it(m_matrix,j); it; ++it)
+        {
+          Index i = it.index();
+          if(i==j)
+          {
+            int k = count[i]++;
+            dest._innerIndexPtr()[k] = i;
+            dest._valuePtr()[k] = it.value();
+          }
+          else if((UpLo==Lower && i>j) || (UpLo==Upper && i<j))
+          {
+            int k = count[j]++;
+            dest._innerIndexPtr()[k] = i;
+            dest._valuePtr()[k] = it.value();
+            k = count[i]++;
+            dest._innerIndexPtr()[k] = j;
+            dest._valuePtr()[k] = internal::conj(it.value());
+          }
+        }
+      }
+    }
+    
+        
 
     // const SparseLLT<PlainObject, UpLo> llt() const;
     // const SparseLDLT<PlainObject, UpLo> ldlt() const;
@@ -100,6 +176,8 @@ template<typename MatrixType, unsigned int UpLo> class SparseSelfAdjointView
   protected:
 
     const typename MatrixType::Nested m_matrix;
+    VectorI m_countPerRow;
+    VectorI m_countPerCol;
 };
 
 /***************************************************************************
@@ -133,7 +211,7 @@ SparseSelfAdjointView<MatrixType,UpLo>::rankUpdate(const SparseMatrixBase<Derive
   if(alpha==Scalar(0))
     m_matrix.const_cast_derived() = tmp.template triangularView<UpLo>();
   else
-    m_matrix.const_cast_derived() /*+*/= alpha * tmp.template triangularView<UpLo>();
+    m_matrix.const_cast_derived() += alpha * tmp.template triangularView<UpLo>();
 
   return *this;
 }
diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp
index 92fff4057..9d79ca740 100644
--- a/test/sparse_basic.cpp
+++ b/test/sparse_basic.cpp
@@ -251,10 +251,21 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re
     VERIFY_IS_APPROX(m2, refM2);
   }
   
+  // test selfadjointView
+  {
+    DenseMatrix refMat2(rows, rows), refMat3(rows, rows);
+    SparseMatrixType m2(rows, rows), m3(rows, rows);
+    initSparse<Scalar>(density, refMat2, m2);
+    refMat3 = refMat2.template selfadjointView<Lower>();
+    m3 = m2.template selfadjointView<Lower>();
+    VERIFY_IS_APPROX(m3, refMat3);
+  }
+  
   // test sparseView
   {
     DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
     SparseMatrixType m2(rows, rows);
+    initSparse<Scalar>(density, refMat2, m2);
     VERIFY_IS_APPROX(m2.eval(), refMat2.sparseView().eval());
   }
 }