From ddb3ac98a20d5f56146a53d485b1899a46b9f912 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Thu, 23 Jul 2009 10:05:38 +0200
Subject: [PATCH] addd matrix * self adjoint high level API

---
 Eigen/src/Core/SelfAdjointView.h | 99 ++++++++++++++++++++------------
 test/product_selfadjoint.cpp     |  7 +++
 2 files changed, 69 insertions(+), 37 deletions(-)

diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h
index 7f5fd7533..c64ebc174 100644
--- a/Eigen/src/Core/SelfAdjointView.h
+++ b/Eigen/src/Core/SelfAdjointView.h
@@ -53,8 +53,9 @@ struct ei_traits<SelfAdjointView<MatrixType, TriangularPart> > : ei_traits<Matri
   };
 };
 
-template<typename Lhs,typename Rhs,bool RhsIsVector=Rhs::IsVectorAtCompileTime>
-struct ei_selfadjoint_matrix_product_returntype;
+template <typename Lhs, int LhsMode, bool LhsIsVector,
+          typename Rhs, int RhsMode, bool RhsIsVector>
+struct ei_selfadjoint_product_returntype;
 
 // FIXME could also be called SelfAdjointWrapper to be consistent with DiagonalWrapper ??
 template<typename MatrixType, unsigned int UpLo> class SelfAdjointView
@@ -99,10 +100,22 @@ template<typename MatrixType, unsigned int UpLo> class SelfAdjointView
 
     /** Efficient self-adjoint matrix times vector/matrix product */
     template<typename OtherDerived>
-    ei_selfadjoint_matrix_product_returntype<SelfAdjointView,OtherDerived>
+    ei_selfadjoint_product_returntype<MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime>
     operator*(const MatrixBase<OtherDerived>& rhs) const
     {
-      return ei_selfadjoint_matrix_product_returntype<SelfAdjointView,OtherDerived>(*this, rhs.derived());
+      return ei_selfadjoint_product_returntype
+              <MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime>
+              (m_matrix, rhs.derived());
+    }
+
+    /** Efficient vector/matrix times self-adjoint matrix product */
+    template<typename OtherDerived> friend
+    ei_selfadjoint_product_returntype<OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false>
+    operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView& rhs)
+    {
+      return ei_selfadjoint_product_returntype
+              <OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false>
+              (lhs.derived(),rhs.m_matrix);
     }
 
     /** Perform a symmetric rank 2 update of the selfadjoint matrix \c *this:
@@ -125,6 +138,14 @@ template<typename MatrixType, unsigned int UpLo> class SelfAdjointView
     const typename MatrixType::Nested m_matrix;
 };
 
+
+// template<typename OtherDerived, typename MatrixType, unsigned int UpLo>
+// ei_selfadjoint_matrix_product_returntype<OtherDerived,SelfAdjointView<MatrixType,UpLo> >
+// operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView<MatrixType,UpLo>& rhs)
+// {
+//   return ei_matrix_selfadjoint_product_returntype<OtherDerived,SelfAdjointView<MatrixType,UpLo> >(lhs.derived(),rhs);
+// }
+
 template<typename Derived1, typename Derived2, int UnrollCount, bool ClearOpposite>
 struct ei_triangular_assignment_selector<Derived1, Derived2, SelfAdjoint, UnrollCount, ClearOpposite>
 {
@@ -163,14 +184,14 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, SelfAdjoint, Dynami
 * Wrapper to ei_product_selfadjoint_vector
 ***************************************************************************/
 
-template<typename Lhs,typename Rhs>
-struct ei_selfadjoint_matrix_product_returntype<Lhs,Rhs,true>
-  : public ReturnByValue<ei_selfadjoint_matrix_product_returntype<Lhs,Rhs,true>,
+template<typename Lhs, int LhsMode, typename Rhs, int RhsMode>
+struct ei_selfadjoint_product_returntype<Lhs,LhsMode,false,Rhs,RhsMode,true>
+  : public ReturnByValue<ei_selfadjoint_product_returntype<Lhs,LhsMode,false,Rhs,RhsMode,true>,
                          Matrix<typename ei_traits<Rhs>::Scalar,
                                 Rhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
 {
   typedef typename ei_cleantype<typename Rhs::Nested>::type RhsNested;
-  ei_selfadjoint_matrix_product_returntype(const Lhs& lhs, const Rhs& rhs)
+  ei_selfadjoint_product_returntype(const Lhs& lhs, const Rhs& rhs)
     : m_lhs(lhs), m_rhs(rhs)
   {}
 
@@ -178,10 +199,10 @@ struct ei_selfadjoint_matrix_product_returntype<Lhs,Rhs,true>
   {
     dst.resize(m_rhs.rows(), m_rhs.cols());
     ei_product_selfadjoint_vector<typename Lhs::Scalar,ei_traits<Lhs>::Flags&RowMajorBit,
-      Lhs::Mode&(UpperTriangularBit|LowerTriangularBit)>
+      LhsMode&(UpperTriangularBit|LowerTriangularBit)>
       (
         m_lhs.rows(), // size
-        m_lhs._expression().data(), // lhs
+        m_lhs.data(), // lhs
         m_lhs.stride(), // lhsStride,
         m_rhs.data(), // rhs
         // int rhsIncr,
@@ -189,7 +210,7 @@ struct ei_selfadjoint_matrix_product_returntype<Lhs,Rhs,true>
       );
   }
 
-  const Lhs m_lhs;
+  const typename Lhs::Nested m_lhs;
   const typename Rhs::Nested m_rhs;
 };
 
@@ -197,25 +218,36 @@ struct ei_selfadjoint_matrix_product_returntype<Lhs,Rhs,true>
 * Wrapper to ei_product_selfadjoint_matrix
 ***************************************************************************/
 
-template<typename Lhs,typename Rhs>
-struct ei_selfadjoint_matrix_product_returntype<Lhs,Rhs,false>
-  : public ReturnByValue<ei_selfadjoint_matrix_product_returntype<Lhs,Rhs,false>,
+template<typename Lhs, int LhsMode, typename Rhs, int RhsMode>
+struct ei_selfadjoint_product_returntype<Lhs,LhsMode,false,Rhs,RhsMode,false>
+  : public ReturnByValue<ei_selfadjoint_product_returntype<Lhs,LhsMode,false,Rhs,RhsMode,false>,
                          Matrix<typename ei_traits<Rhs>::Scalar,
-                                Rhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
+                                Lhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
 {
-  ei_selfadjoint_matrix_product_returntype(const Lhs& lhs, const Rhs& rhs)
+  ei_selfadjoint_product_returntype(const Lhs& lhs, const Rhs& rhs)
     : m_lhs(lhs), m_rhs(rhs)
   {}
 
   typedef typename Lhs::Scalar Scalar;
+
+  typedef typename Lhs::Nested LhsNested;
+  typedef typename ei_cleantype<LhsNested>::type _LhsNested;
+  typedef ei_blas_traits<_LhsNested> LhsBlasTraits;
+  typedef typename LhsBlasTraits::DirectLinearAccessType ActualLhsType;
+  typedef typename ei_cleantype<ActualLhsType>::type _ActualLhsType;
+  
   typedef typename Rhs::Nested RhsNested;
   typedef typename ei_cleantype<RhsNested>::type _RhsNested;
+  typedef ei_blas_traits<_RhsNested> RhsBlasTraits;
+  typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhsType;
+  typedef typename ei_cleantype<ActualRhsType>::type _ActualRhsType;
 
-  typedef typename ei_traits<Lhs>::ExpressionType LhsExpr;
-  typedef typename LhsExpr::Nested LhsNested;
-  typedef typename ei_cleantype<LhsNested>::type _LhsNested;
-
-  enum { UpLo = ei_traits<Lhs>::Mode&(UpperTriangularBit|LowerTriangularBit) };
+  enum {
+    LhsUpLo = LhsMode&(UpperTriangularBit|LowerTriangularBit),
+    LhsIsSelfAdjoint = (LhsMode&SelfAdjointBit)==SelfAdjointBit,
+    RhsUpLo = RhsMode&(UpperTriangularBit|LowerTriangularBit),
+    RhsIsSelfAdjoint = (RhsMode&SelfAdjointBit)==SelfAdjointBit
+  };
 
   template<typename Dest> inline void _addTo(Dest& dst) const
   { evalTo(dst,1); }
@@ -231,26 +263,19 @@ struct ei_selfadjoint_matrix_product_returntype<Lhs,Rhs,false>
 
   template<typename Dest> void evalTo(Dest& dst, Scalar alpha) const
   {
-    typedef ei_blas_traits<_LhsNested> LhsBlasTraits;
-    typedef ei_blas_traits<_RhsNested> RhsBlasTraits;
-
-    typedef typename LhsBlasTraits::DirectLinearAccessType ActualLhsType;
-    typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhsType;
-
-    typedef typename ei_cleantype<ActualLhsType>::type _ActualLhsType;
-    typedef typename ei_cleantype<ActualRhsType>::type _ActualRhsType;
-
-    const ActualLhsType lhs = LhsBlasTraits::extract(m_lhs._expression());
+    const ActualLhsType lhs = LhsBlasTraits::extract(m_lhs);
     const ActualRhsType rhs = RhsBlasTraits::extract(m_rhs);
 
-    Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(m_lhs._expression())
+    Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(m_lhs)
                                * RhsBlasTraits::extractScalarFactor(m_rhs);
 
     ei_product_selfadjoint_matrix<Scalar,
-      EIGEN_LOGICAL_XOR(UpLo==UpperTriangular,
-      ei_traits<Lhs>::Flags &RowMajorBit) ? RowMajor : ColMajor, true,
-      NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(UpLo==UpperTriangular,bool(LhsBlasTraits::NeedToConjugate)),
-      ei_traits<Rhs>::Flags &RowMajorBit  ? RowMajor : ColMajor, false, bool(RhsBlasTraits::NeedToConjugate),
+      EIGEN_LOGICAL_XOR(LhsUpLo==UpperTriangular,
+                        ei_traits<Lhs>::Flags &RowMajorBit) ? RowMajor : ColMajor, LhsIsSelfAdjoint,
+      NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(LhsUpLo==UpperTriangular,bool(LhsBlasTraits::NeedToConjugate)),
+      EIGEN_LOGICAL_XOR(RhsUpLo==UpperTriangular,
+                        ei_traits<Rhs>::Flags &RowMajorBit) ? RowMajor : ColMajor, RhsIsSelfAdjoint,
+      NumTraits<Scalar>::IsComplex && EIGEN_LOGICAL_XOR(RhsUpLo==UpperTriangular,bool(RhsBlasTraits::NeedToConjugate)),
       ei_traits<Dest>::Flags&RowMajorBit  ? RowMajor : ColMajor>
       ::run(
         lhs.rows(), rhs.cols(),   // sizes
@@ -261,7 +286,7 @@ struct ei_selfadjoint_matrix_product_returntype<Lhs,Rhs,false>
       );
   }
 
-  const Lhs m_lhs;
+  const LhsNested m_lhs;
   const RhsNested m_rhs;
 };
 
diff --git a/test/product_selfadjoint.cpp b/test/product_selfadjoint.cpp
index 814d542e4..44bafad93 100644
--- a/test/product_selfadjoint.cpp
+++ b/test/product_selfadjoint.cpp
@@ -138,6 +138,13 @@ template<typename MatrixType> void symm(const MatrixType& m)
   m2 = m1.template triangularView<UpperTriangular>();
   VERIFY_IS_APPROX(rhs32 = (s1*m2.adjoint()).template selfadjointView<LowerTriangular>() * (s2*rhs3).conjugate(),
                    rhs33 = (s1*m1.adjoint()) * (s2*rhs3).conjugate());
+
+  // test matrix * selfadjoint
+  m2 = m1.template triangularView<LowerTriangular>();
+  VERIFY_IS_APPROX(rhs22 = (rhs2) * (m2).template selfadjointView<LowerTriangular>(),
+                   rhs23 = (rhs2) * (m1));
+  VERIFY_IS_APPROX(rhs22 = (s2*rhs2) * (s1*m2).template selfadjointView<LowerTriangular>(),
+                   rhs23 = (s2*rhs2) * (s1*m1));
 }
 void test_product_selfadjoint()
 {