add a generic mechanism to copy a special matrix to a dense matrix so that

we don't need to add other specialization of MatrixBase::operator=, Matrix::=,
and Matrix::Matrix(...)

6 files changed, 53 insertions, 112 deletions

diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
index bad360579..fccd53459 100644
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h
@@ -225,47 +225,6 @@ template<typename MatrixType> struct LLT_Traits<MatrixType,UpperTriangular>
 template<typename MatrixType, int _UpLo>
 void LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
 {
-//   assert(a.rows()==a.cols());
-//   const int size = a.rows();
-//   m_matrix.resize(size, size);
-//   // The biggest overall is the point of reference to which further diagonals
-//   // are compared; if any diagonal is negligible compared
-//   // to the largest overall, the algorithm bails.  This cutoff is suggested
-//   // in "Analysis of the Cholesky Decomposition of a Semi-definite Matrix" by
-//   // Nicholas J. Higham. Also see "Accuracy and Stability of Numerical
-//   // Algorithms" page 217, also by Higham.
-//   const RealScalar cutoff = machine_epsilon<Scalar>() * size * a.diagonal().cwise().abs().maxCoeff();
-//   RealScalar x;
-//   x = ei_real(a.coeff(0,0));
-//   m_matrix.coeffRef(0,0) = ei_sqrt(x);
-//   if(size==1)
-//   {
-//     m_isInitialized = true;
-//     return;
-//   }
-//   m_matrix.col(0).end(size-1) = a.row(0).end(size-1).adjoint() / ei_real(m_matrix.coeff(0,0));
-//   for (int j = 1; j < size; ++j)
-//   {
-//     x = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).squaredNorm();
-//     if (ei_abs(x) < cutoff) continue;
-// 
-//     m_matrix.coeffRef(j,j) = x = ei_sqrt(x);
-// 
-//     int endSize = size-j-1;
-//     if (endSize>0) {
-//       // Note that when all matrix columns have good alignment, then the following
-//       // product is guaranteed to be optimal with respect to alignment.
-//       m_matrix.col(j).end(endSize) =
-//         (m_matrix.block(j+1, 0, endSize, j) * m_matrix.row(j).start(j).adjoint()).lazy();
-// 
-//       // FIXME could use a.col instead of a.row
-//       m_matrix.col(j).end(endSize) = (a.row(j).end(endSize).adjoint()
-//         - m_matrix.col(j).end(endSize) ) / x;
-//     }
-//   }
-// 
-//   m_isInitialized = true;
-
   assert(a.rows()==a.cols());
   const int size = a.rows();
   m_matrix.resize(size, size);
diff --git a/Eigen/src/Core/DiagonalMatrix.h b/Eigen/src/Core/DiagonalMatrix.h
index a092b7794..5fc80c92b 100644
--- a/Eigen/src/Core/DiagonalMatrix.h
+++ b/Eigen/src/Core/DiagonalMatrix.h
@@ -50,6 +50,8 @@ class DiagonalBase : public MultiplierBase<Derived>
     #endif // not EIGEN_PARSED_BY_DOXYGEN
     
     DenseMatrixType toDenseMatrix() const { return derived(); }
+    template<typename DenseDerived>
+    void evalToDense(MatrixBase<DenseDerived> &other) const;
 
     inline const DiagonalVectorType& diagonal() const { return derived().diagonal(); }
     inline DiagonalVectorType& diagonal() { return derived().diagonal(); }
@@ -63,12 +65,11 @@ class DiagonalBase : public MultiplierBase<Derived>
 };
 
 template<typename Derived>
-template<typename DiagonalDerived>
-Derived& MatrixBase<Derived>::operator=(const DiagonalBase<DiagonalDerived> &other)
+template<typename DenseDerived>
+void DiagonalBase<Derived>::evalToDense(MatrixBase<DenseDerived> &other) const
 {
-  setZero();
-  diagonal() = other.diagonal();
-  return derived();
+  other.setZero();
+  other.diagonal() = diagonal();
 }
 
 /** \class DiagonalMatrix
diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h
index b7a8ebf21..18c8f969a 100644
--- a/Eigen/src/Core/Matrix.h
+++ b/Eigen/src/Core/Matrix.h
@@ -439,34 +439,20 @@ class Matrix
     { other.evalTo(*this); }
     /** Destructor */
     inline ~Matrix() {}
-    
-
-    template<typename DiagonalDerived>
-    EIGEN_STRONG_INLINE Matrix& operator=(const DiagonalBase<DiagonalDerived> &other)
-    {
-      resize(other.diagonal().size(), other.diagonal().size());
-      Base::operator=(other);
-      return *this;
-    }
 
-    template<typename DiagonalDerived>
-    EIGEN_STRONG_INLINE Matrix(const DiagonalBase<DiagonalDerived> &other)
-      : m_storage(other.diagonal().size() * other.diagonal().size(), other.diagonal().size(), other.diagonal().size())
-    {
-      *this = other;
-    }
-
-    template<typename TriangularDerived>
-    EIGEN_STRONG_INLINE Matrix& operator=(const TriangularBase<TriangularDerived> &other)
+    /** \sa MatrixBase::operator=(const AnyMatrixBase<OtherDerived>&) */
+    template<typename OtherDerived>
+    EIGEN_STRONG_INLINE Matrix& operator=(const AnyMatrixBase<OtherDerived> &other)
     {
-      resize(other.rows(), other.cols());
+      resize(other.derived().rows(), other.derived().cols());
       Base::operator=(other.derived());
       return *this;
     }
 
-    template<typename TriangularDerived>
-    EIGEN_STRONG_INLINE Matrix(const TriangularBase<TriangularDerived> &other)
-      : m_storage(other.rows() * other.cols(), other.rows(), other.cols())
+    /** \sa MatrixBase::operator=(const AnyMatrixBase<OtherDerived>&) */
+    template<typename OtherDerived>
+    EIGEN_STRONG_INLINE Matrix(const AnyMatrixBase<OtherDerived> &other)
+      : m_storage(other.derived().rows() * other.derived().cols(), other.derived().rows(), other.derived().cols())
     {
       *this = other;
     }
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index 82417a144..7aedbd3c7 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -284,6 +284,17 @@ template<typename Derived> class MatrixBase
       */
     Derived& operator=(const MatrixBase& other);
 
+    /** Copies the generic expression \a other into *this. \returns a reference to *this.
+      * The expression must provide a (templated) evalToDense(Derived& dst) const function
+      * which does the actual job. In practice, this allows any user to write its own
+      * special matrix without having to modify MatrixBase */
+    template<typename OtherDerived>
+    Derived& operator=(const AnyMatrixBase<OtherDerived> &other)
+    { other.derived().evalToDense(derived()); return derived(); }
+
+    template<typename OtherDerived,typename OtherEvalType>
+    Derived& operator=(const ReturnByValue<OtherDerived,OtherEvalType>& func);
+
 #ifndef EIGEN_PARSED_BY_DOXYGEN
     /** Copies \a other into *this without evaluating other. \returns a reference to *this. */
     template<typename OtherDerived>
@@ -297,10 +308,6 @@ template<typename Derived> class MatrixBase
     template<typename OtherDerived>
     Derived& lazyAssign(const Flagged<OtherDerived, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit>& other)
     { return lazyAssign(other._expression()); }
-
-    /** Overloaded for fast triangular part to dense matrix evaluation */
-    template<typename TriangularDerived>
-    Derived& lazyAssign(const TriangularBase<TriangularDerived> &other);
 #endif // not EIGEN_PARSED_BY_DOXYGEN
 
     CommaInitializer<Derived> operator<< (const Scalar& s);
@@ -403,12 +410,6 @@ template<typename Derived> class MatrixBase
     const DiagonalProduct<Derived, DiagonalDerived, DiagonalOnTheRight>
     operator*(const DiagonalBase<DiagonalDerived> &diagonal) const;
 
-    template<typename DiagonalDerived>
-    Derived& operator=(const DiagonalBase<DiagonalDerived> &other);
-
-    template<typename TriangularDerived>
-    Derived& operator=(const TriangularBase<TriangularDerived> &other);
-
     template<typename OtherDerived>
     typename ei_plain_matrix_type_column_major<OtherDerived>::type
     solveTriangular(const MatrixBase<OtherDerived>& other) const;
@@ -738,9 +739,6 @@ template<typename Derived> class MatrixBase
     // dense = dense * sparse
     template<typename Derived1, typename Derived2>
     Derived& lazyAssign(const SparseProduct<Derived1,Derived2,DenseTimeSparseProduct>& product);
-
-    template<typename OtherDerived,typename OtherEvalType>
-    Derived& operator=(const ReturnByValue<OtherDerived,OtherEvalType>& func);
     
     #ifdef EIGEN_MATRIXBASE_PLUGIN
     #include EIGEN_MATRIXBASE_PLUGIN
diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h
index 8131ef323..81621dbcc 100644
--- a/Eigen/src/Core/TriangularMatrix.h
+++ b/Eigen/src/Core/TriangularMatrix.h
@@ -90,6 +90,11 @@ template<typename Derived> class TriangularBase : public MultiplierBase<Derived>
     inline Derived& derived() { return *static_cast<Derived*>(this); }
     #endif // not EIGEN_PARSED_BY_DOXYGEN
 
+    template<typename DenseDerived>
+    void evalToDense(MatrixBase<DenseDerived> &other) const;
+    template<typename DenseDerived>
+    void evalToDenseLazy(MatrixBase<DenseDerived> &other) const;
+
   protected:
 
     void check_coordinates(int row, int col)
@@ -516,36 +521,34 @@ void TriangularView<MatrixType, Mode>::lazyAssign(const TriangularBase<OtherDeri
 /** Assigns a triangular or selfadjoint matrix to a dense matrix.
   * If the matrix is triangular, the opposite part is set to zero. */
 template<typename Derived>
-template<typename TriangularDerived>
-Derived& MatrixBase<Derived>::operator=(const TriangularBase<TriangularDerived> &other)
+template<typename DenseDerived>
+void TriangularBase<Derived>::evalToDense(MatrixBase<DenseDerived> &other) const
 {
-  if(ei_traits<TriangularDerived>::Flags & EvalBeforeAssigningBit)
+  if(ei_traits<Derived>::Flags & EvalBeforeAssigningBit)
   {
-    typename TriangularDerived::PlainMatrixType other_evaluated(other.rows(), other.cols());
-    other_evaluated.lazyAssign(other);
-    this->swap(other_evaluated);
+    typename Derived::PlainMatrixType other_evaluated(rows(), cols());
+    evalToDenseLazy(other_evaluated);
+    other.derived().swap(other_evaluated);
   }
   else
-    lazyAssign(other.derived());
-  return derived();
+    evalToDenseLazy(other.derived());
 }
 
 /** Assigns a triangular or selfadjoint matrix to a dense matrix.
   * If the matrix is triangular, the opposite part is set to zero. */
 template<typename Derived>
-template<typename TriangularDerived>
-Derived& MatrixBase<Derived>::lazyAssign(const TriangularBase<TriangularDerived> &other)
+template<typename DenseDerived>
+void TriangularBase<Derived>::evalToDenseLazy(MatrixBase<DenseDerived> &other) const
 {
-  const bool unroll =   Derived::SizeAtCompileTime * TriangularDerived::CoeffReadCost / 2
+  const bool unroll =   DenseDerived::SizeAtCompileTime * Derived::CoeffReadCost / 2
                      <= EIGEN_UNROLLING_LIMIT;
   ei_assert(this->rows() == other.rows() && this->cols() == other.cols());
 
   ei_part_assignment_impl
-    <Derived, typename ei_traits<TriangularDerived>::ExpressionType, TriangularDerived::Mode,
-    unroll ? int(Derived::SizeAtCompileTime) : Dynamic,
+    <DenseDerived, typename ei_traits<Derived>::ExpressionType, Derived::Mode,
+    unroll ? int(DenseDerived::SizeAtCompileTime) : Dynamic,
     true // clear the opposite triangular part
-    >::run(this->const_cast_derived(), other.derived()._expression());
-  return derived();
+    >::run(other.derived(), derived()._expression());
 }
 
 /** \deprecated use MatrixBase::triangularView() */
diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h
index 8c53757b0..82527f7ef 100644
--- a/Eigen/src/Sparse/SparseMatrixBase.h
+++ b/Eigen/src/Sparse/SparseMatrixBase.h
@@ -451,25 +451,19 @@ template<typename Derived> class SparseMatrixBase
 //     Derived& setRandom();
 //     Derived& setIdentity();
 
-      Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime> toDense() const
+      template<typename DenseDerived>
+      void evalToDense(MatrixBase<DenseDerived>& dst)
       {
-        Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime> res(rows(),cols());
-        res.setZero();
+        dst.resize(rows(),cols());
+        dst.setZero();
         for (int j=0; j<outerSize(); ++j)
-        {
           for (typename Derived::InnerIterator i(derived(),j); i; ++i)
-          {
-            if(IsRowMajor)
-              std::cerr << i.row() << "," << i.col() << " == " << j << "," << i.index() << "\n";
-            else
-              std::cerr << i.row() << "," << i.col() << " == " << i.index() << "," << j << "\n";
-//             if(IsRowMajor)
-              res.coeffRef(i.row(),i.col()) = i.value();
-//             else
-//               res.coeffRef(i.index(),j) = i.value();
-          }
-        }
-        return res;
+            res.coeffRef(i.row(),i.col()) = i.value();
+      }
+
+      Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime> toDense() const
+      {
+        return derived();
       }
 
     template<typename OtherDerived>