Attachment 916 Details for Bug 1658 – Add conjugateIf<cond>() members

This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen

[patch] Add conjugateIf<cond>() members

add_conjugate_if.diff (text/plain), 10.82 KB, created by Gael Guennebaud on 2019-01-15 09:58:26 UTC

(hide)

Description:

Filename:

MIME Type:

Creator: Gael Guennebaud

Created: 2019-01-15 09:58:26 UTC

Size: 10.82 KB

patch

obsolete

>diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h
>--- a/Eigen/src/Core/SelfAdjointView.h
>+++ b/Eigen/src/Core/SelfAdjointView.h
>@@ -56,16 +56,17 @@ template<typename _MatrixType, unsigned 
>     typedef typename internal::traits<SelfAdjointView>::MatrixTypeNested MatrixTypeNested;
>     typedef typename internal::traits<SelfAdjointView>::MatrixTypeNestedCleaned MatrixTypeNestedCleaned;
>     typedef MatrixTypeNestedCleaned NestedExpression;
> 
>     /** \brief The type of coefficients in this matrix */
>     typedef typename internal::traits<SelfAdjointView>::Scalar Scalar; 
>     typedef typename MatrixType::StorageIndex StorageIndex;
>     typedef typename internal::remove_all<typename MatrixType::ConjugateReturnType>::type MatrixConjugateReturnType;
>+    typedef SelfAdjointView<typename internal::add_const<MatrixType>::type, UpLo> ConstSelfAdjointView;
> 
>     enum {
>       Mode = internal::traits<SelfAdjointView>::Mode,
>       Flags = internal::traits<SelfAdjointView>::Flags,
>       TransposeMode = ((Mode & Upper) ? Lower : 0) | ((Mode & Lower) ? Upper : 0)
>     };
>     typedef typename MatrixType::PlainObject PlainObject;
> 
>@@ -192,16 +193,28 @@ template<typename _MatrixType, unsigned 
>     }
> 
>     typedef SelfAdjointView<const MatrixConjugateReturnType,UpLo> ConjugateReturnType;
>     /** \sa MatrixBase::conjugate() const */
>     EIGEN_DEVICE_FUNC
>     inline const ConjugateReturnType conjugate() const
>     { return ConjugateReturnType(m_matrix.conjugate()); }
> 
>+    /** \returns an expression of the complex conjugate of \c *this if Cond==true,
>+     *           returns \c *this otherwise.
>+     */
>+    template<bool Cond>
>+    EIGEN_DEVICE_FUNC
>+    inline typename internal::conditional<Cond,ConjugateReturnType,ConstSelfAdjointView>::type
>+    conjugateIf() const
>+    {
>+      typedef typename internal::conditional<Cond,ConjugateReturnType,ConstSelfAdjointView>::type ReturnType;
>+      return ReturnType(m_matrix.template conjugateIf<Cond>());
>+    }
>+
>     typedef SelfAdjointView<const typename MatrixType::AdjointReturnType,TransposeMode> AdjointReturnType;
>     /** \sa MatrixBase::adjoint() const */
>     EIGEN_DEVICE_FUNC
>     inline const AdjointReturnType adjoint() const
>     { return AdjointReturnType(m_matrix.adjoint()); }
> 
>     typedef SelfAdjointView<typename MatrixType::TransposeReturnType,TransposeMode> TransposeReturnType;
>      /** \sa MatrixBase::transpose() */
>diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h
>--- a/Eigen/src/Core/TriangularMatrix.h
>+++ b/Eigen/src/Core/TriangularMatrix.h
>@@ -193,16 +193,17 @@ template<typename _MatrixType, unsigned 
>     typedef typename internal::traits<TriangularView>::Scalar Scalar;
>     typedef _MatrixType MatrixType;
> 
>   protected:
>     typedef typename internal::traits<TriangularView>::MatrixTypeNested MatrixTypeNested;
>     typedef typename internal::traits<TriangularView>::MatrixTypeNestedNonRef MatrixTypeNestedNonRef;
> 
>     typedef typename internal::remove_all<typename MatrixType::ConjugateReturnType>::type MatrixConjugateReturnType;
>+    typedef TriangularView<typename internal::add_const<MatrixType>::type, _Mode> ConstTriangularView;
>     
>   public:
> 
>     typedef typename internal::traits<TriangularView>::StorageKind StorageKind;
>     typedef typename internal::traits<TriangularView>::MatrixTypeNestedCleaned NestedExpression;
> 
>     enum {
>       Mode = _Mode,
>@@ -238,16 +239,28 @@ template<typename _MatrixType, unsigned 
>     NestedExpression& nestedExpression() { return m_matrix; }
>     
>     typedef TriangularView<const MatrixConjugateReturnType,Mode> ConjugateReturnType;
>     /** \sa MatrixBase::conjugate() const */
>     EIGEN_DEVICE_FUNC
>     inline const ConjugateReturnType conjugate() const
>     { return ConjugateReturnType(m_matrix.conjugate()); }
> 
>+    /** \returns an expression of the complex conjugate of \c *this if Cond==true,
>+     *           returns \c *this otherwise.
>+     */
>+    template<bool Cond>
>+    EIGEN_DEVICE_FUNC
>+    inline typename internal::conditional<Cond,ConjugateReturnType,ConstTriangularView>::type
>+    conjugateIf() const
>+    {
>+      typedef typename internal::conditional<Cond,ConjugateReturnType,ConstTriangularView>::type ReturnType;
>+      return ReturnType(m_matrix.template conjugateIf<Cond>());
>+    }
>+
>     typedef TriangularView<const typename MatrixType::AdjointReturnType,TransposeMode> AdjointReturnType;
>     /** \sa MatrixBase::adjoint() const */
>     EIGEN_DEVICE_FUNC
>     inline const AdjointReturnType adjoint() const
>     { return AdjointReturnType(m_matrix.adjoint()); }
> 
>     typedef TriangularView<typename MatrixType::TransposeReturnType,TransposeMode> TransposeReturnType;
>      /** \sa MatrixBase::transpose() */
>diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h
>--- a/Eigen/src/LU/PartialPivLU.h
>+++ b/Eigen/src/LU/PartialPivLU.h
>@@ -241,36 +241,31 @@ template<typename _MatrixType> class Par
> 
>       // Step 3
>       m_lu.template triangularView<Upper>().solveInPlace(dst);
>     }
> 
>     template<bool Conjugate, typename RhsType, typename DstType>
>     EIGEN_DEVICE_FUNC
>     void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const {
>-     /* The decomposition PA = LU can be rewritten as A = P^{-1} L U.
>+     /* The decomposition PA = LU can be rewritten as A^T = U^T L^T P.
>       * So we proceed as follows:
>-      * Step 1: compute c = Pb.
>-      * Step 2: replace c by the solution x to Lx = c.
>-      * Step 3: replace c by the solution x to Ux = c.
>+      * Step 1: compute c as the solution to L^T c = b
>+      * Step 2: replace c by the solution x to U^T x = c.
>+      * Step 3: update  c = P^-1 c.
>       */
> 
>       eigen_assert(rhs.rows() == m_lu.cols());
> 
>-      if (Conjugate) {
>-        // Step 1
>-        dst = m_lu.template triangularView<Upper>().adjoint().solve(rhs);
>-        // Step 2
>-        m_lu.template triangularView<UnitLower>().adjoint().solveInPlace(dst);
>-      } else {
>-        // Step 1
>-        dst = m_lu.template triangularView<Upper>().transpose().solve(rhs);
>-        // Step 2
>-        m_lu.template triangularView<UnitLower>().transpose().solveInPlace(dst);
>-      }
>+      // Step 1
>+      dst = m_lu.template triangularView<Upper>().transpose()
>+                .template conjugateIf<Conjugate>().solve(rhs);
>+      // Step 2
>+      m_lu.template triangularView<UnitLower>().transpose()
>+          .template conjugateIf<Conjugate>().solveInPlace(dst);
>       // Step 3
>       dst = permutationP().transpose() * dst;
>     }
>     #endif
> 
>   protected:
> 
>     static void check_template_parameters()
>diff --git a/Eigen/src/plugins/CommonCwiseUnaryOps.h b/Eigen/src/plugins/CommonCwiseUnaryOps.h
>--- a/Eigen/src/plugins/CommonCwiseUnaryOps.h
>+++ b/Eigen/src/plugins/CommonCwiseUnaryOps.h
>@@ -71,16 +71,30 @@ EIGEN_DOC_UNARY_ADDONS(conjugate,complex
> /// \sa <a href="group__CoeffwiseMathFunctions.html#cwisetable_conj">Math functions</a>, MatrixBase::adjoint()
> EIGEN_DEVICE_FUNC
> inline ConjugateReturnType
> conjugate() const
> {
>   return ConjugateReturnType(derived());
> }
> 
>+/// \returns an expression of the complex conjugate of \c *this if Cond==true, returns derived() otherwise.
>+///
>+EIGEN_DOC_UNARY_ADDONS(conjugate,complex conjugate)
>+///
>+/// \sa conjugate()
>+template<bool Cond>
>+EIGEN_DEVICE_FUNC
>+inline typename internal::conditional<Cond,ConjugateReturnType,const Derived&>::type
>+conjugateIf() const
>+{
>+  typedef typename internal::conditional<Cond,ConjugateReturnType,const Derived&>::type ReturnType;
>+  return ReturnType(derived());
>+}
>+
> /// \returns a read-only expression of the real part of \c *this.
> ///
> EIGEN_DOC_UNARY_ADDONS(real,real part function)
> ///
> /// \sa imag()
> EIGEN_DEVICE_FUNC
> inline RealReturnType
> real() const { return RealReturnType(derived()); }
>diff --git a/test/adjoint.cpp b/test/adjoint.cpp
>--- a/test/adjoint.cpp
>+++ b/test/adjoint.cpp
>@@ -138,16 +138,19 @@ template<typename MatrixType> void adjoi
>   m3.transposeInPlace();
>   VERIFY_IS_APPROX(m3,m1.conjugate());
> 
>   // check mixed dot product
>   typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealVectorType;
>   RealVectorType rv1 = RealVectorType::Random(rows);
>   VERIFY_IS_APPROX(v1.dot(rv1.template cast<Scalar>()), v1.dot(rv1));
>   VERIFY_IS_APPROX(rv1.template cast<Scalar>().dot(v1), rv1.dot(v1));
>+
>+  VERIFY( is_same_type(m1,m1.template conjugateIf<false>()) );
>+  VERIFY( is_same_type(m1.conjugate(),m1.template conjugateIf<true>()) );
> }
> 
> template<int>
> void adjoint_extra()
> {
>   MatrixXcf a(10,10), b(10,10);
>   VERIFY_RAISES_ASSERT(a = a.transpose());
>   VERIFY_RAISES_ASSERT(a = a.transpose() + b);
>diff --git a/test/triangular.cpp b/test/triangular.cpp
>--- a/test/triangular.cpp
>+++ b/test/triangular.cpp
>@@ -124,16 +124,32 @@ template<typename MatrixType> void trian
>   m1up = m1.template triangularView<Upper>();
>   VERIFY_IS_APPROX(m1.template selfadjointView<Upper>().template triangularView<Upper>().toDenseMatrix(), m1up);
>   VERIFY_IS_APPROX(m1up.template selfadjointView<Upper>().template triangularView<Upper>().toDenseMatrix(), m1up);
>   VERIFY_IS_APPROX(m1.template selfadjointView<Upper>().template triangularView<Lower>().toDenseMatrix(), m1up.adjoint());
>   VERIFY_IS_APPROX(m1up.template selfadjointView<Upper>().template triangularView<Lower>().toDenseMatrix(), m1up.adjoint());
> 
>   VERIFY_IS_APPROX(m1.template selfadjointView<Upper>().diagonal(), m1.diagonal());
> 
>+  m3.setRandom();
>+  const MatrixType& m3c(m3);
>+  VERIFY( is_same_type(m3c.template triangularView<Lower>(),m3.template triangularView<Lower>().template conjugateIf<false>()) );
>+  VERIFY( is_same_type(m3c.template triangularView<Lower>().conjugate(),m3.template triangularView<Lower>().template conjugateIf<true>()) );
>+  VERIFY_IS_APPROX(m3.template triangularView<Lower>().template conjugateIf<true>().toDenseMatrix(),
>+                   m3.conjugate().template triangularView<Lower>().toDenseMatrix());
>+  VERIFY_IS_APPROX(m3.template triangularView<Lower>().template conjugateIf<false>().toDenseMatrix(),
>+                   m3.template triangularView<Lower>().toDenseMatrix());
>+
>+  VERIFY( is_same_type(m3c.template selfadjointView<Lower>(),m3.template selfadjointView<Lower>().template conjugateIf<false>()) );
>+  VERIFY( is_same_type(m3c.template selfadjointView<Lower>().conjugate(),m3.template selfadjointView<Lower>().template conjugateIf<true>()) );
>+  VERIFY_IS_APPROX(m3.template selfadjointView<Lower>().template conjugateIf<true>().toDenseMatrix(),
>+                   m3.conjugate().template selfadjointView<Lower>().toDenseMatrix());
>+  VERIFY_IS_APPROX(m3.template selfadjointView<Lower>().template conjugateIf<false>().toDenseMatrix(),
>+                   m3.template selfadjointView<Lower>().toDenseMatrix());
>+
> }
> 
> 
> template<typename MatrixType> void triangular_rect(const MatrixType& m)
> {
>   typedef typename MatrixType::Scalar Scalar;
>   typedef typename NumTraits<Scalar>::Real RealScalar;
>   enum { Rows =  MatrixType::RowsAtCompileTime, Cols =  MatrixType::ColsAtCompileTime };

Actions: View | Diff

Attachments on bug 1658: 916