added optimized matrix times diagonal matrix product via Diagonal flag shortcut.

8 files changed, 56 insertions, 17 deletions

diff --git a/Eigen/src/Array/CwiseOperators.h b/Eigen/src/Array/CwiseOperators.h
index 82a18d11c..f9f9e1267 100644
--- a/Eigen/src/Array/CwiseOperators.h
+++ b/Eigen/src/Array/CwiseOperators.h
@@ -76,6 +76,14 @@ MatrixBase<Derived>::cwisePow(const Scalar& exponent) const
     (derived(), ei_scalar_pow_op<Scalar>(exponent));
 }
 
+/** \returns an expression of the coefficient-wise reciprocal of *this. */
+template<typename Derived>
+inline const CwiseUnaryOp<ei_scalar_inverse_op<typename ei_traits<Derived>::Scalar>, Derived>
+MatrixBase<Derived>::cwiseInverse() const
+{
+  return derived();
+}
+
 // -- binary operators --
 
 /** \returns an expression of the coefficient-wise \< operator of *this and \a other
diff --git a/Eigen/src/Array/Functors.h b/Eigen/src/Array/Functors.h
index 805f2dcad..3169106d8 100644
--- a/Eigen/src/Array/Functors.h
+++ b/Eigen/src/Array/Functors.h
@@ -100,6 +100,19 @@ template<typename Scalar>
 struct ei_functor_traits<ei_scalar_pow_op<Scalar> >
 { enum { Cost = 5 * NumTraits<Scalar>::MulCost, IsVectorizable = false }; };
 
+/** \internal
+  * \brief Template functor to compute the reciprocal of a scalar
+  *
+  * \sa class CwiseUnaryOp, MatrixBase::cwiseInverse
+  */
+template<typename Scalar>
+struct ei_scalar_inverse_op {
+  inline Scalar operator() (const Scalar& a) const { return Scalar(1)/a; }
+};
+template<typename Scalar>
+struct ei_functor_traits<ei_scalar_inverse_op<Scalar> >
+{ enum { Cost = NumTraits<Scalar>::MulCost, IsVectorizable = false }; };
+
 // default ei_functor_traits for STL functors:
 
 template<typename T>
diff --git a/Eigen/src/Core/DiagonalMatrix.h b/Eigen/src/Core/DiagonalMatrix.h
index 0581c669c..fa343e7d7 100644
--- a/Eigen/src/Core/DiagonalMatrix.h
+++ b/Eigen/src/Core/DiagonalMatrix.h
@@ -49,7 +49,7 @@ struct ei_traits<DiagonalMatrix<CoeffsVectorType> >
     ColsAtCompileTime = CoeffsVectorType::SizeAtCompileTime,
     MaxRowsAtCompileTime = CoeffsVectorType::MaxSizeAtCompileTime,
     MaxColsAtCompileTime = CoeffsVectorType::MaxSizeAtCompileTime,
-    Flags = _CoeffsVectorTypeTypeNested::Flags & HereditaryBits,
+    Flags = _CoeffsVectorTypeTypeNested::Flags & HereditaryBits | Diagonal,
     CoeffReadCost = _CoeffsVectorTypeTypeNested::CoeffReadCost
   };
 };
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index fc1d524bd..ab35b94a6 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -544,6 +544,7 @@ template<typename Derived> class MatrixBase : public ArrayBase<Derived>
     const CwiseUnaryOp<ei_scalar_sin_op<typename ei_traits<Derived>::Scalar>, Derived> cwiseSin() const;
     const CwiseUnaryOp<ei_scalar_pow_op<typename ei_traits<Derived>::Scalar>, Derived>
     cwisePow(const Scalar& exponent) const;
+    const CwiseUnaryOp<ei_scalar_inverse_op<typename ei_traits<Derived>::Scalar>, Derived> cwiseInverse() const;
 
     template<typename OtherDerived>
     const CwiseBinaryOp<std::less<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h
index 1def64777..7643ac610 100644
--- a/Eigen/src/Core/Product.h
+++ b/Eigen/src/Core/Product.h
@@ -147,6 +147,7 @@ template<typename Lhs, typename Rhs> struct ei_product_eval_mode
   enum{ value =  Lhs::MaxRowsAtCompileTime >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
               && Rhs::MaxColsAtCompileTime >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
               && Lhs::MaxColsAtCompileTime >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
+              && Rhs::Flags&Diagonal!=Diagonal
               ? CacheFriendlyProduct : NormalProduct };
 };
 
@@ -259,25 +260,40 @@ template<typename Lhs, typename Rhs, int EvalMode> class Product : ei_no_assignm
 
     const Scalar _coeff(int row, int col) const
     {
-      Scalar res;
-      const bool unroll = CoeffReadCost <= EIGEN_UNROLLING_LIMIT;
-      ei_product_impl<Lhs::ColsAtCompileTime-1,
-                          unroll ? Lhs::ColsAtCompileTime : Dynamic,
-                          _LhsNested, _RhsNested>
-        ::run(row, col, m_lhs, m_rhs, res);
-      return res;
+      if ((Rhs::Flags&Diagonal)==Diagonal)
+      {
+        return  m_lhs.coeff(row, col) * m_rhs.coeff(col, col);
+      }
+      else
+      {
+        Scalar res;
+        const bool unroll = CoeffReadCost <= EIGEN_UNROLLING_LIMIT;
+        ei_product_impl<Lhs::ColsAtCompileTime-1,
+                            unroll ? Lhs::ColsAtCompileTime : Dynamic,
+                            _LhsNested, _RhsNested>
+          ::run(row, col, m_lhs, m_rhs, res);
+        return res;
+      }
     }
 
     template<int LoadMode>
     const PacketScalar _packetCoeff(int row, int col) const
     {
-      const bool unroll = CoeffReadCost <= EIGEN_UNROLLING_LIMIT;
-      PacketScalar res;
-      ei_packet_product_impl<Flags&RowMajorBit ? true : false, Lhs::ColsAtCompileTime-1,
-                          unroll ? Lhs::ColsAtCompileTime : Dynamic,
-                          _LhsNested, _RhsNested, PacketScalar>
-        ::run(row, col, m_lhs, m_rhs, res);
-      return res;
+      if ((Rhs::Flags&Diagonal)==Diagonal)
+      {
+        assert(_LhsNested::Flags&RowMajorBit==0);
+        return ei_pmul(m_lhs.template packetCoeff<LoadMode>(row, col), ei_pset1(m_rhs.coeff(col, col)));
+      }
+      else
+      {
+        const bool unroll = CoeffReadCost <= EIGEN_UNROLLING_LIMIT;
+        PacketScalar res;
+        ei_packet_product_impl<Flags&RowMajorBit ? true : false, Lhs::ColsAtCompileTime-1,
+                            unroll ? Lhs::ColsAtCompileTime : Dynamic,
+                            _LhsNested, _RhsNested, PacketScalar>
+          ::run(row, col, m_lhs, m_rhs, res);
+        return res;
+      }
     }
 
     template<typename Lhs_, typename Rhs_, int EvalMode_, typename DestDerived_, bool DirectAccess_>
diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h
index 909921f35..e17563c9a 100644
--- a/Eigen/src/Core/util/Constants.h
+++ b/Eigen/src/Core/util/Constants.h
@@ -66,6 +66,7 @@ const unsigned int SelfAdjoint = SelfAdjointBit;
 // additional possible values for the Mode parameter of extract()
 const unsigned int UnitUpper = UpperTriangularBit | UnitDiagBit;
 const unsigned int UnitLower = LowerTriangularBit | UnitDiagBit;
+const unsigned int Diagonal = Upper | Lower;
 
 
 
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index cc759ddea..9fa3718d3 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -67,6 +67,7 @@ template<typename Scalar> struct ei_scalar_log_op;
 template<typename Scalar> struct ei_scalar_cos_op;
 template<typename Scalar> struct ei_scalar_sin_op;
 template<typename Scalar> struct ei_scalar_pow_op;
+template<typename Scalar> struct ei_scalar_inverse_op;
 template<typename Scalar, typename NewType> struct ei_scalar_cast_op;
 template<typename Scalar, bool IsVectorizable> struct ei_scalar_multiple_op;
 template<typename Scalar> struct ei_scalar_quotient1_op;
diff --git a/Eigen/src/Core/util/Meta.h b/Eigen/src/Core/util/Meta.h
index f1939f59e..1ce7f3cfe 100644
--- a/Eigen/src/Core/util/Meta.h
+++ b/Eigen/src/Core/util/Meta.h
@@ -212,8 +212,7 @@ template<typename T, int n=1> struct ei_nested
 
 template<unsigned int Flags> struct ei_are_flags_consistent
 {
-  enum { ret = !( (Flags&UnitDiagBit && Flags&ZeroDiagBit)
-                  || (Flags&UpperTriangularBit && Flags&LowerTriangularBit) )
+  enum { ret = !( (Flags&UnitDiagBit && Flags&ZeroDiagBit) )
   };
 };