Refactor triangular assignment

Eigen/src/Core/AssignEvaluator.h

@@ -140,6 +140,7 @@ public:
 template<typename Kernel, int Index, int Stop>
 struct copy_using_evaluator_DefaultTraversal_CompleteUnrolling
 {
+  // FIXME: this is not very clean, perhaps this information should be provided by the kernel?
   typedef typename Kernel::DstEvaluatorType DstEvaluatorType;
   typedef typename DstEvaluatorType::XprType DstXprType;
   
@@ -204,6 +205,7 @@ struct copy_using_evaluator_LinearTraversal_CompleteUnrolling<Kernel, Stop, Stop
 template<typename Kernel, int Index, int Stop>
 struct copy_using_evaluator_innervec_CompleteUnrolling
 {
+  // FIXME: this is not very clean, perhaps this information should be provided by the kernel?
   typedef typename Kernel::DstEvaluatorType DstEvaluatorType;
   typedef typename DstEvaluatorType::XprType DstXprType;
   
@@ -233,8 +235,7 @@ struct copy_using_evaluator_innervec_InnerUnrolling
   static EIGEN_STRONG_INLINE void run(Kernel &kernel, int outer)
   {
     kernel.template assignPacketByOuterInner<Aligned, Aligned>(outer, Index);
-    typedef typename Kernel::DstEvaluatorType::XprType DstXprType;
+    enum { NextIndex = Index + packet_traits<typename Kernel::Scalar>::size };
-    enum { NextIndex = Index + packet_traits<typename DstXprType::Scalar>::size };
     copy_using_evaluator_innervec_InnerUnrolling<Kernel, NextIndex, Stop>::run(kernel, outer);
   }
 };
@@ -542,12 +543,6 @@ public:
     m_functor.assignCoeff(m_dst.coeffRef(row,col), m_src.coeff(row,col));
   }
   
-  /// This overload by-passes the source expression, i.e., dst(row,col) ?= value
-  void assignCoeff(Index row, Index col, const Scalar &value)
-  {
-    m_functor.assignCoeff(m_dst.coeffRef(row,col), value);
-  }
-  
   /// \sa assignCoeff(Index,Index)
   void assignCoeff(Index index)
   {

Eigen/src/Core/SelfAdjointView.h

@@ -244,6 +244,8 @@ template<typename _MatrixType, unsigned int UpLo> class SelfAdjointView
 
 namespace internal {
 
+#ifndef EIGEN_TEST_EVALUATORS
+ 
 template<typename Derived1, typename Derived2, int UnrollCount, bool ClearOpposite>
 struct triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Upper), UnrollCount, ClearOpposite>
 {
@@ -336,6 +338,8 @@ struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Lower, Dyn
   }
 };
 
+#endif // EIGEN_TEST_EVALUATORS
+
 #ifdef EIGEN_ENABLE_EVALUATORS
 // TODO currently a selfadjoint expression has the form SelfAdjointView<.,.>
 //      in the future selfadjoint-ness should be defined by the expression traits
@@ -348,6 +352,7 @@ struct evaluator_traits<SelfAdjointView<MatrixType,Mode> >
   
   static const int AssumeAliasing = 0;
 };
+
 #endif // EIGEN_ENABLE_EVALUATORS
 
 } // end namespace internal

Eigen/src/Core/Swap.h

@@ -148,6 +148,7 @@ template<typename DstEvaluatorTypeT, typename SrcEvaluatorTypeT>
 class generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, swap_assign_op<typename DstEvaluatorTypeT::Scalar>, Specialized>
  : public generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, swap_assign_op<typename DstEvaluatorTypeT::Scalar>, BuiltIn>
 {
+protected:
   typedef generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, swap_assign_op<typename DstEvaluatorTypeT::Scalar>, BuiltIn> Base;
   typedef typename DstEvaluatorTypeT::PacketScalar PacketScalar;
   using Base::m_dst;

Eigen/src/Core/TriangularMatrix.h

@@ -469,6 +469,8 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView
         return m_matrix.diagonal().prod();
     }
     
+#ifndef EIGEN_TEST_EVALUATORS
+    
     // TODO simplify the following:
     template<typename ProductDerived, typename Lhs, typename Rhs>
     EIGEN_DEVICE_FUNC
@@ -515,6 +517,8 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView
       return assignProduct(other,-other.alpha());
     }
 
+#endif // EIGEN_TEST_EVALUATORS
+
   protected:
     
     template<typename ProductDerived, typename Lhs, typename Rhs>
@@ -530,6 +534,8 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView
 
 namespace internal {
 
+#ifndef EIGEN_TEST_EVALUATORS
+  
 template<typename Derived1, typename Derived2, unsigned int Mode, int UnrollCount, bool ClearOpposite>
 struct triangular_assignment_selector
 {
@@ -694,8 +700,52 @@ struct triangular_assignment_selector<Derived1, Derived2, UnitLower, Dynamic, Cl
   }
 };
 
+#endif // EIGEN_TEST_EVALUATORS
+
 } // end namespace internal
 
+#ifdef EIGEN_TEST_EVALUATORS
+
+// FIXME should we keep that possibility
+template<typename MatrixType, unsigned int Mode>
+template<typename OtherDerived>
+inline TriangularView<MatrixType, Mode>&
+TriangularView<MatrixType, Mode>::operator=(const MatrixBase<OtherDerived>& other)
+{
+  internal::call_assignment(*this, other.template triangularView<Mode>(), internal::assign_op<Scalar>());
+  return *this;
+}
+
+// FIXME should we keep that possibility
+template<typename MatrixType, unsigned int Mode>
+template<typename OtherDerived>
+void TriangularView<MatrixType, Mode>::lazyAssign(const MatrixBase<OtherDerived>& other)
+{
+  internal::call_assignment(this->noalias(), other.template triangularView<Mode>());
+}
+
+
+
+template<typename MatrixType, unsigned int Mode>
+template<typename OtherDerived>
+inline TriangularView<MatrixType, Mode>&
+TriangularView<MatrixType, Mode>::operator=(const TriangularBase<OtherDerived>& other)
+{
+  eigen_assert(Mode == int(OtherDerived::Mode));
+  internal::call_assignment(*this, other.derived());
+  return *this;
+}
+
+template<typename MatrixType, unsigned int Mode>
+template<typename OtherDerived>
+void TriangularView<MatrixType, Mode>::lazyAssign(const TriangularBase<OtherDerived>& other)
+{
+  eigen_assert(Mode == int(OtherDerived::Mode));
+  internal::call_assignment(this->noalias(), other.derived());
+}
+
+#else
+
 // FIXME should we keep that possibility
 template<typename MatrixType, unsigned int Mode>
 template<typename OtherDerived>
@@ -770,6 +820,8 @@ void TriangularView<MatrixType, Mode>::lazyAssign(const TriangularBase<OtherDeri
     >::run(m_matrix.const_cast_derived(), other.derived().nestedExpression());
 }
 
+#endif // EIGEN_TEST_EVALUATORS
+
 /***************************************************************************
 * Implementation of TriangularBase methods
 ***************************************************************************/
@@ -790,6 +842,8 @@ void TriangularBase<Derived>::evalTo(MatrixBase<DenseDerived> &other) const
     evalToLazy(other.derived());
 }
 
+#ifndef EIGEN_TEST_EVALUATORS
+
 /** 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>
@@ -811,6 +865,8 @@ void TriangularBase<Derived>::evalToLazy(MatrixBase<DenseDerived> &other) const
     >::run(other.derived(), derived().nestedExpression());
 }
 
+#endif // EIGEN_TEST_EVALUATORS
+
 /***************************************************************************
 * Implementation of TriangularView methods
 ***************************************************************************/
@@ -962,15 +1018,15 @@ struct evaluator_traits<TriangularView<MatrixType,Mode> >
 
 template<typename MatrixType, unsigned int Mode, typename Kind>
 struct evaluator<TriangularView<MatrixType,Mode>, Kind, typename MatrixType::Scalar>
- : evaluator<MatrixType>
+ : evaluator<typename internal::remove_all<MatrixType>::type>
 {
   typedef TriangularView<MatrixType,Mode> XprType;
-  typedef evaluator<MatrixType> Base;
+  typedef evaluator<typename internal::remove_all<MatrixType>::type> Base;
   typedef evaluator type;
   evaluator(const XprType &xpr) : Base(xpr.nestedExpression()) {}
 };
 
-// Additional assignement kinds:
+// Additional assignment kinds:
 struct Triangular2Triangular    {};
 struct Triangular2Dense         {};
 struct Dense2Triangular         {};
@@ -979,6 +1035,59 @@ struct Dense2Triangular         {};
 template<typename Kernel, unsigned int Mode, int UnrollCount, bool ClearOpposite> struct triangular_assignment_loop;
 
  
+// Specialization of the dense assignment kernel for triangular matrices.
+// The main additions are:
+//  - An overload of assignCoeff(i, j, val) that by-passes the source expression. Needed to set the opposite triangular part to 0.
+//  - When calling assignCoeff(i,j...), we guarantee that i!=j
+//  - New assignDiagonalCoeff(i), and assignDiagonalCoeff(i, val) for assigning coefficients on the diagonal.
+template<int Mode, int ClearOpposite, typename DstEvaluatorTypeT, typename SrcEvaluatorTypeT, typename Functor, int Version = Specialized>
+class triangular_dense_assignment_kernel : public generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version>
+{
+protected:
+  typedef generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version> Base;
+  typedef typename Base::DstXprType DstXprType;
+  typedef typename Base::SrcXprType SrcXprType;
+  using Base::m_dst;
+  using Base::m_src;
+  using Base::m_functor;
+public:
+  
+  typedef typename Base::DstEvaluatorType DstEvaluatorType;
+  typedef typename Base::SrcEvaluatorType SrcEvaluatorType;
+  typedef typename Base::Scalar Scalar;
+  typedef typename Base::Index Index;
+  typedef typename Base::AssignmentTraits AssignmentTraits;
+  
+  
+  triangular_dense_assignment_kernel(DstEvaluatorType &dst, const SrcEvaluatorType &src, const Functor &func, DstXprType& dstExpr)
+    : Base(dst, src, func, dstExpr)
+  {}
+  
+#ifdef EIGEN_INTERNAL_DEBUGGING
+  void assignCoeff(Index row, Index col)
+  {
+    eigen_internal_assert(row!=col);
+    Base::assignCoeff(row,col);
+  }
+#else
+  using Base::assignCoeff;
+#endif
+  
+  void assignDiagonalCoeff(Index id)
+  {
+         if((Mode&UnitDiag) && ClearOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(1));
+    else if((Mode&ZeroDiag) && ClearOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(0));
+    else if((Mode&(UnitDiag|ZeroDiag))==0)    Base::assignCoeff(id,id);
+  }
+  
+  void assignOppositeCoeff(Index row, Index col)
+  { 
+    eigen_internal_assert(row!=col);
+    if(ClearOpposite)
+      m_functor.assignCoeff(m_dst.coeffRef(row,col), Scalar(0));
+  }
+};
+
 template<int Mode, bool ClearOpposite, typename DstXprType, typename SrcXprType, typename Functor>
 void call_triangular_assignment_loop(const DstXprType& dst, const SrcXprType& src, const Functor &func)
 {
@@ -994,13 +1103,13 @@ void call_triangular_assignment_loop(const DstXprType& dst, const SrcXprType& sr
   DstEvaluatorType dstEvaluator(dst);
   SrcEvaluatorType srcEvaluator(src);
     
-  typedef generic_dense_assignment_kernel<DstEvaluatorType,SrcEvaluatorType,Functor> Kernel;
+  typedef triangular_dense_assignment_kernel<Mode,ClearOpposite,DstEvaluatorType,SrcEvaluatorType,Functor> Kernel;
   Kernel kernel(dstEvaluator, srcEvaluator, func, dst.const_cast_derived());
   
   enum {
       unroll = DstXprType::SizeAtCompileTime != Dynamic
-            && internal::traits<SrcXprType>::CoeffReadCost != Dynamic
+//             && internal::traits<SrcXprType>::CoeffReadCost != Dynamic
-            && DstXprType::SizeAtCompileTime * internal::traits<SrcXprType>::CoeffReadCost / 2 <= EIGEN_UNROLLING_LIMIT
+//             && DstXprType::SizeAtCompileTime * internal::traits<SrcXprType>::CoeffReadCost / 2 <= EIGEN_UNROLLING_LIMIT
     };
   
   triangular_assignment_loop<Kernel, Mode, unroll ? int(DstXprType::SizeAtCompileTime) : Dynamic, ClearOpposite>::run(kernel);
@@ -1050,9 +1159,13 @@ struct Assignment<DstXprType, SrcXprType, Functor, Dense2Triangular, Scalar>
 template<typename Kernel, unsigned int Mode, int UnrollCount, bool ClearOpposite>
 struct triangular_assignment_loop
 {
+  // FIXME: this is not very clean, perhaps this information should be provided by the kernel?
+  typedef typename Kernel::DstEvaluatorType DstEvaluatorType;
+  typedef typename DstEvaluatorType::XprType DstXprType;
+  
   enum {
-    col = (UnrollCount-1) / Kernel::DstEvaluatorType::RowsAtCompileTime,
+    col = (UnrollCount-1) / DstXprType::RowsAtCompileTime,
-    row = (UnrollCount-1) % Kernel::DstEvaluatorType::RowsAtCompileTime
+    row = (UnrollCount-1) % DstXprType::RowsAtCompileTime
   };
   
   typedef typename Kernel::Scalar Scalar;
@@ -1062,23 +1175,12 @@ struct triangular_assignment_loop
   {
     triangular_assignment_loop<Kernel, Mode, UnrollCount-1, ClearOpposite>::run(kernel);
     
-    // TODO should be a static assert
+    if(row==col)
-    eigen_assert(    Mode == Upper              || Mode == Lower
+      kernel.assignDiagonalCoeff(row);
-                  || Mode == StrictlyUpper      || Mode == StrictlyLower
+    else if( ((Mode&Lower) && row>col) || ((Mode&Upper) && row<col) )
-                  || Mode == UnitUpper          || Mode == UnitLower);
+      kernel.assignCoeff(row,col);
-    
+    else
-    if((Mode == Upper           && row <= col)
+      kernel.assignOppositeCoeff(row,col);
-    || (Mode == Lower           && row >= col)
-    || (Mode == StrictlyUpper   && row <  col)
-    || (Mode == StrictlyLower   && row >  col)
-    || (Mode == UnitUpper       && row <  col)
-    || (Mode == UnitLower       && row >  col))
-      kernel.assignCoeff(row, col);
-    else if(ClearOpposite)
-    {
-      if((Mode&UnitDiag) && row==col) kernel.assignCoeff(row, col, Scalar(1));
-      else                            kernel.assignCoeff(row, col, Scalar(0));
-    }
   }
 };
 
@@ -1090,94 +1192,14 @@ struct triangular_assignment_loop<Kernel, Mode, 0, ClearOpposite>
   static inline void run(Kernel &) {}
 };
 
-// TODO: merge the following 6 variants into a single one (or max two),
-//       and perhaps write them using sub-expressions
 
-// TODO: expreiment with a recursive assignment procedure splitting the current
+
+// TODO: experiment with a recursive assignment procedure splitting the current
 //       triangular part into one rectangular and two triangular parts.
 
-template<typename Kernel, bool ClearOpposite>
-struct triangular_assignment_loop<Kernel, Upper, Dynamic, ClearOpposite>
-{
-  typedef typename Kernel::Index Index;
-  typedef typename Kernel::Scalar Scalar;
-  EIGEN_DEVICE_FUNC
-  static inline void run(Kernel &kernel)
-  {
-    for(Index j = 0; j < kernel.cols(); ++j)
-    {
-      Index maxi = (std::min)(j, kernel.rows()-1);
-      for(Index i = 0; i <= maxi; ++i)
-        kernel.assignCoeff(i, j);
-      if (ClearOpposite)
-        for(Index i = maxi+1; i < kernel.rows(); ++i)
-          kernel.assignCoeff(i, j, Scalar(0));
-    }
-  }
-};
 
-template<typename Kernel, bool ClearOpposite>
+template<typename Kernel, unsigned int Mode, bool ClearOpposite>
-struct triangular_assignment_loop<Kernel, Lower, Dynamic, ClearOpposite>
+struct triangular_assignment_loop<Kernel, Mode, Dynamic, ClearOpposite>
-{
-  typedef typename Kernel::Index Index;
-  typedef typename Kernel::Scalar Scalar;
-  EIGEN_DEVICE_FUNC
-  static inline void run(Kernel &kernel)
-  {
-    for(Index j = 0; j < kernel.cols(); ++j)
-    {
-      for(Index i = j; i < kernel.rows(); ++i)
-        kernel.assignCoeff(i, j);
-      Index maxi = (std::min)(j, kernel.rows());
-      if (ClearOpposite)
-        for(Index i = 0; i < maxi; ++i)
-          kernel.assignCoeff(i, j, Scalar(0));
-    }
-  }
-};
-
-template<typename Kernel, bool ClearOpposite>
-struct triangular_assignment_loop<Kernel, StrictlyUpper, Dynamic, ClearOpposite>
-{
-  typedef typename Kernel::Index Index;
-  typedef typename Kernel::Scalar Scalar;
-  EIGEN_DEVICE_FUNC
-  static inline void run(Kernel &kernel)
-  {
-    for(Index j = 0; j < kernel.cols(); ++j)
-    {
-      Index maxi = (std::min)(j, kernel.rows());
-      for(Index i = 0; i < maxi; ++i)
-        kernel.assignCoeff(i, j);
-      if (ClearOpposite)
-        for(Index i = maxi; i < kernel.rows(); ++i)
-          kernel.assignCoeff(i, j, Scalar(0));
-    }
-  }
-};
-
-template<typename Kernel, bool ClearOpposite>
-struct triangular_assignment_loop<Kernel, StrictlyLower, Dynamic, ClearOpposite>
-{
-  typedef typename Kernel::Index Index;
-  typedef typename Kernel::Scalar Scalar;
-  EIGEN_DEVICE_FUNC
-  static inline void run(Kernel &kernel)
-  {
-    for(Index j = 0; j < kernel.cols(); ++j)
-    {
-      for(Index i = j+1; i < kernel.rows(); ++i)
-        kernel.assignCoeff(i, j);
-      Index maxi = (std::min)(j, kernel.rows()-1);
-      if (ClearOpposite)
-        for(Index i = 0; i <= maxi; ++i)
-          kernel.assignCoeff(i, j, Scalar(0));
-    }
-  }
-};
-
-template<typename Kernel, bool ClearOpposite>
-struct triangular_assignment_loop<Kernel, UnitUpper, Dynamic, ClearOpposite>
 {
   typedef typename Kernel::Index Index;
   typedef typename Kernel::Scalar Scalar;
@@ -1188,50 +1210,41 @@ struct triangular_assignment_loop<Kernel, UnitUpper, Dynamic, ClearOpposite>
     {
       Index maxi = (std::min)(j, kernel.rows());
       Index i = 0;
-      for(; i < maxi; ++i)
+      if (((Mode&Lower) && ClearOpposite) || (Mode&Upper))
-        kernel.assignCoeff(i, j);
-      
-      if(i<kernel.rows()) // then i==j
-        kernel.assignCoeff(i++, j, Scalar(1));
-      
-      if (ClearOpposite)
-      {
-        for(; i < kernel.rows(); ++i)
-          kernel.assignCoeff(i, j, Scalar(0));
-      }
-    }
-  }
-};
-template<typename Kernel, bool ClearOpposite>
-struct triangular_assignment_loop<Kernel, UnitLower, Dynamic, ClearOpposite>
-{
-  typedef typename Kernel::Index Index;
-  typedef typename Kernel::Scalar Scalar;
-  EIGEN_DEVICE_FUNC
-  static inline void run(Kernel &kernel)
-  {
-    for(Index j = 0; j < kernel.cols(); ++j)
-    {
-      Index maxi = (std::min)(j, kernel.rows());
-      Index i = 0;
-      if (ClearOpposite)
       {
         for(; i < maxi; ++i)
-          kernel.assignCoeff(i, j, Scalar(0));
+          if(Mode&Upper) kernel.assignCoeff(i, j);
+          else           kernel.assignOppositeCoeff(i, j);
       }
+      else
+        i = maxi;
       
       if(i<kernel.rows()) // then i==j
-        kernel.assignCoeff(i++, j, Scalar(1));
+        kernel.assignDiagonalCoeff(i++);
       
+      if (((Mode&Upper) && ClearOpposite) || (Mode&Lower))
+      {
         for(; i < kernel.rows(); ++i)
-        kernel.assignCoeff(i, j);
+          if(Mode&Lower) kernel.assignCoeff(i, j);
+          else           kernel.assignOppositeCoeff(i, j);
+      }
     }
   }
 };
 
 } // end namespace internal
 
-#endif
+/** 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 DenseDerived>
+void TriangularBase<Derived>::evalToLazy(MatrixBase<DenseDerived> &other) const
+{
+  other.derived().resize(this->rows(), this->cols());
+  internal::call_triangular_assignment_loop<Derived::Mode,true /* ClearOpposite */>(other.derived(), derived().nestedExpression());
+}
+
+#endif // EIGEN_ENABLE_EVALUATORS
 
 } // end namespace Eigen