From 5b78780deffd4a5728cfc6ae5371f56e038e2ac0 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Tue, 18 Feb 2014 17:43:16 +0100
Subject: [PATCH] Add evaluator shortcut for triangular ?= product

---
 Eigen/src/Core/TriangularMatrix.h             | 70 ++++++++++++++++++-
 .../products/GeneralMatrixMatrixTriangular.h  | 12 +++-
 test/product_notemporary.cpp                  |  3 +-
 3 files changed, 79 insertions(+), 6 deletions(-)

diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h
index c5f137ab8..ffa2faeaa 100644
--- a/Eigen/src/Core/TriangularMatrix.h
+++ b/Eigen/src/Core/TriangularMatrix.h
@@ -233,6 +233,24 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView
     EIGEN_DEVICE_FUNC
     inline Index innerStride() const { return m_matrix.innerStride(); }
 
+#ifdef EIGEN_TEST_EVALUATORS
+
+    /** \sa MatrixBase::operator+=() */    
+    template<typename Other>
+    EIGEN_DEVICE_FUNC
+    TriangularView&  operator+=(const DenseBase<Other>& other) {
+      internal::call_assignment_no_alias(*this, other.derived(), internal::add_assign_op<Scalar>());
+      return *this;
+    }
+    /** \sa MatrixBase::operator-=() */
+    template<typename Other>
+    EIGEN_DEVICE_FUNC
+    TriangularView&  operator-=(const DenseBase<Other>& other) {
+      internal::call_assignment_no_alias(*this, other.derived(), internal::sub_assign_op<Scalar>());
+      return *this;
+    }
+    
+#else
     /** \sa MatrixBase::operator+=() */    
     template<typename Other>
     EIGEN_DEVICE_FUNC
@@ -241,6 +259,7 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView
     template<typename Other>
     EIGEN_DEVICE_FUNC
     TriangularView&  operator-=(const DenseBase<Other>& other) { return *this = m_matrix - other.derived(); }
+#endif
     /** \sa MatrixBase::operator*=() */
     EIGEN_DEVICE_FUNC
     TriangularView&  operator*=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = m_matrix * other; }
@@ -527,12 +546,18 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView
 
 #endif // EIGEN_TEST_EVALUATORS
 
-  protected:
+#ifdef EIGEN_TEST_EVALUATORS
+    template<typename ProductType>
+    EIGEN_DEVICE_FUNC
+    EIGEN_STRONG_INLINE TriangularView& _assignProduct(const ProductType& prod, const Scalar& alpha);
     
+  protected:
+#else
+  protected:
     template<typename ProductDerived, typename Lhs, typename Rhs>
     EIGEN_DEVICE_FUNC
     EIGEN_STRONG_INLINE TriangularView& assignProduct(const ProductBase<ProductDerived, Lhs,Rhs>& prod, const Scalar& alpha);
-
+#endif
     MatrixTypeNested m_matrix;
 };
 
@@ -720,7 +745,7 @@ 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>());
+  internal::call_assignment_no_alias(*this, other.derived(), internal::assign_op<Scalar>());
   return *this;
 }
 
@@ -1253,6 +1278,45 @@ void TriangularBase<Derived>::evalToLazy(MatrixBase<DenseDerived> &other) const
   other.derived().resize(this->rows(), this->cols());
   internal::call_triangular_assignment_loop<Derived::Mode,(Derived::Mode&SelfAdjoint)==0 /* SetOpposite */>(other.derived(), derived().nestedExpression());
 }
+
+namespace internal {
+  
+// Triangular = Product
+template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
+struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::assign_op<Scalar>, Dense2Triangular, Scalar>
+{
+  typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
+  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &)
+  {
+    dst.setZero();
+    dst._assignProduct(src, 1);
+  }
+};
+
+// Triangular += Product
+template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
+struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::add_assign_op<Scalar>, Dense2Triangular, Scalar>
+{
+  typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
+  static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op<Scalar> &)
+  {
+    dst._assignProduct(src, 1);
+  }
+};
+
+// Triangular -= Product
+template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar>
+struct Assignment<DstXprType, Product<Lhs,Rhs,DefaultProduct>, internal::sub_assign_op<Scalar>, Dense2Triangular, Scalar>
+{
+  typedef Product<Lhs,Rhs,DefaultProduct> SrcXprType;
+  static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op<Scalar> &)
+  {
+    dst._assignProduct(src, -1);
+  }
+};
+
+
+} // end namespace internal
 #endif
 
 #endif // EIGEN_ENABLE_EVALUATORS
diff --git a/Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h b/Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h
index 87bd36bc3..46be46f85 100644
--- a/Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h
+++ b/Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h
@@ -264,6 +264,16 @@ struct general_product_to_triangular_selector<MatrixType,ProductType,UpLo,false>
   }
 };
 
+#ifdef EIGEN_TEST_EVALUATORS
+template<typename MatrixType, unsigned int UpLo>
+template<typename ProductType>
+TriangularView<MatrixType,UpLo>& TriangularView<MatrixType,UpLo>::_assignProduct(const ProductType& prod, const Scalar& alpha)
+{
+  general_product_to_triangular_selector<MatrixType, ProductType, UpLo, internal::traits<ProductType>::InnerSize==1>::run(m_matrix.const_cast_derived(), prod, alpha);
+  
+  return *this;
+}
+#else
 template<typename MatrixType, unsigned int UpLo>
 template<typename ProductDerived, typename _Lhs, typename _Rhs>
 TriangularView<MatrixType,UpLo>& TriangularView<MatrixType,UpLo>::assignProduct(const ProductBase<ProductDerived, _Lhs,_Rhs>& prod, const Scalar& alpha)
@@ -272,7 +282,7 @@ TriangularView<MatrixType,UpLo>& TriangularView<MatrixType,UpLo>::assignProduct(
   
   return *this;
 }
-
+#endif
 } // end namespace Eigen
 
 #endif // EIGEN_GENERAL_MATRIX_MATRIX_TRIANGULAR_H
diff --git a/test/product_notemporary.cpp b/test/product_notemporary.cpp
index 258d238e2..deb494c96 100644
--- a/test/product_notemporary.cpp
+++ b/test/product_notemporary.cpp
@@ -114,8 +114,7 @@ template<typename MatrixType> void product_notemporary(const MatrixType& m)
   VERIFY_EVALUATION_COUNT( Scalar tmp = 0; tmp += Scalar(RealScalar(1)) /  (m3.transpose() * m3).diagonal().array().abs().sum(), 0 );
 
   // Zero temporaries for ... CoeffBasedProductMode
-  // - does not work with GCC because of the <..>, we'ld need variadic macros ...
-  //VERIFY_EVALUATION_COUNT( m3.col(0).head<5>() * m3.col(0).transpose() + m3.col(0).head<5>() * m3.col(0).transpose(), 0 );
+  VERIFY_EVALUATION_COUNT( m3.col(0).template head<5>() * m3.col(0).transpose() + m3.col(0).template head<5>() * m3.col(0).transpose(), 0 );
 
   // Check matrix * vectors
   VERIFY_EVALUATION_COUNT( cvres.noalias() = m1 * cv1, 0 );