From 7f32109c11b9cbc3cedc72e59683bf5839d35d75 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 17 Jan 2019 11:33:43 +0100 Subject: [PATCH] Add conjugateIf members to DesneBase, TriangularView, SelfadjointView, and make PartialPivLU use it. --- Eigen/src/Core/SelfAdjointView.h | 13 +++++++++++++ Eigen/src/Core/TriangularMatrix.h | 13 +++++++++++++ Eigen/src/LU/PartialPivLU.h | 25 ++++++++++--------------- Eigen/src/plugins/CommonCwiseUnaryOps.h | 14 ++++++++++++++ test/adjoint.cpp | 3 +++ test/triangular.cpp | 16 ++++++++++++++++ 6 files changed, 69 insertions(+), 15 deletions(-) diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h index 2cf3fa1ef..2173799d9 100644 --- a/Eigen/src/Core/SelfAdjointView.h +++ b/Eigen/src/Core/SelfAdjointView.h @@ -61,6 +61,7 @@ template class SelfAdjointView typedef typename internal::traits::Scalar Scalar; typedef typename MatrixType::StorageIndex StorageIndex; typedef typename internal::remove_all::type MatrixConjugateReturnType; + typedef SelfAdjointView::type, UpLo> ConstSelfAdjointView; enum { Mode = internal::traits::Mode, @@ -197,6 +198,18 @@ template class SelfAdjointView 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 + EIGEN_DEVICE_FUNC + inline typename internal::conditional::type + conjugateIf() const + { + typedef typename internal::conditional::type ReturnType; + return ReturnType(m_matrix.template conjugateIf()); + } + typedef SelfAdjointView AdjointReturnType; /** \sa MatrixBase::adjoint() const */ EIGEN_DEVICE_FUNC diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index 521de6160..cf3532f06 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -198,6 +198,7 @@ template class TriangularView typedef typename internal::traits::MatrixTypeNestedNonRef MatrixTypeNestedNonRef; typedef typename internal::remove_all::type MatrixConjugateReturnType; + typedef TriangularView::type, _Mode> ConstTriangularView; public: @@ -243,6 +244,18 @@ template class TriangularView 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 + EIGEN_DEVICE_FUNC + inline typename internal::conditional::type + conjugateIf() const + { + typedef typename internal::conditional::type ReturnType; + return ReturnType(m_matrix.template conjugateIf()); + } + typedef TriangularView AdjointReturnType; /** \sa MatrixBase::adjoint() const */ EIGEN_DEVICE_FUNC diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index ecc0e748f..ff4be360e 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -246,26 +246,21 @@ template class PartialPivLU template 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().adjoint().solve(rhs); - // Step 2 - m_lu.template triangularView().adjoint().solveInPlace(dst); - } else { - // Step 1 - dst = m_lu.template triangularView().transpose().solve(rhs); - // Step 2 - m_lu.template triangularView().transpose().solveInPlace(dst); - } + // Step 1 + dst = m_lu.template triangularView().transpose() + .template conjugateIf().solve(rhs); + // Step 2 + m_lu.template triangularView().transpose() + .template conjugateIf().solveInPlace(dst); // Step 3 dst = permutationP().transpose() * dst; } diff --git a/Eigen/src/plugins/CommonCwiseUnaryOps.h b/Eigen/src/plugins/CommonCwiseUnaryOps.h index 89f4faaac..5418dc415 100644 --- a/Eigen/src/plugins/CommonCwiseUnaryOps.h +++ b/Eigen/src/plugins/CommonCwiseUnaryOps.h @@ -76,6 +76,20 @@ 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 +EIGEN_DEVICE_FUNC +inline typename internal::conditional::type +conjugateIf() const +{ + typedef typename internal::conditional::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) diff --git a/test/adjoint.cpp b/test/adjoint.cpp index e2bfa6d7d..4c4f98bb9 100644 --- a/test/adjoint.cpp +++ b/test/adjoint.cpp @@ -143,6 +143,9 @@ template void adjoint(const MatrixType& m) RealVectorType rv1 = RealVectorType::Random(rows); VERIFY_IS_APPROX(v1.dot(rv1.template cast()), v1.dot(rv1)); VERIFY_IS_APPROX(rv1.template cast().dot(v1), rv1.dot(v1)); + + VERIFY( is_same_type(m1,m1.template conjugateIf()) ); + VERIFY( is_same_type(m1.conjugate(),m1.template conjugateIf()) ); } template diff --git a/test/triangular.cpp b/test/triangular.cpp index 99ef1dcda..0fca5e3b9 100644 --- a/test/triangular.cpp +++ b/test/triangular.cpp @@ -129,6 +129,22 @@ template void triangular_square(const MatrixType& m) VERIFY_IS_APPROX(m1.template selfadjointView().diagonal(), m1.diagonal()); + m3.setRandom(); + const MatrixType& m3c(m3); + VERIFY( is_same_type(m3c.template triangularView(),m3.template triangularView().template conjugateIf()) ); + VERIFY( is_same_type(m3c.template triangularView().conjugate(),m3.template triangularView().template conjugateIf()) ); + VERIFY_IS_APPROX(m3.template triangularView().template conjugateIf().toDenseMatrix(), + m3.conjugate().template triangularView().toDenseMatrix()); + VERIFY_IS_APPROX(m3.template triangularView().template conjugateIf().toDenseMatrix(), + m3.template triangularView().toDenseMatrix()); + + VERIFY( is_same_type(m3c.template selfadjointView(),m3.template selfadjointView().template conjugateIf()) ); + VERIFY( is_same_type(m3c.template selfadjointView().conjugate(),m3.template selfadjointView().template conjugateIf()) ); + VERIFY_IS_APPROX(m3.template selfadjointView().template conjugateIf().toDenseMatrix(), + m3.conjugate().template selfadjointView().toDenseMatrix()); + VERIFY_IS_APPROX(m3.template selfadjointView().template conjugateIf().toDenseMatrix(), + m3.template selfadjointView().toDenseMatrix()); + }