mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-24 14:45:14 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add conjugateIf<bool> members to DesneBase, TriangularView, SelfadjointView, and make PartialPivLU use it.

Browse Source
This commit is contained in:
Gael Guennebaud 2019-01-17 11:33:43 +01:00
parent 7b35c26b1c
commit 7f32109c11
6 changed files with 69 additions and 15 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
Eigen/src/Core/SelfAdjointView.h
View File

@ -61,6 +61,7 @@ template<typename _MatrixType, unsigned int UpLo> class SelfAdjointView
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,
@ -197,6 +198,18 @@ template<typename _MatrixType, unsigned int UpLo> 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<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

13
Eigen/src/Core/TriangularMatrix.h
View File

@ -198,6 +198,7 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView
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:
@ -243,6 +244,18 @@ template<typename _MatrixType, unsigned int _Mode> 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<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

25
Eigen/src/LU/PartialPivLU.h
View File

@ -246,26 +246,21 @@ template<typename _MatrixType> class PartialPivLU
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;
}

14
Eigen/src/plugins/CommonCwiseUnaryOps.h
View File

@ -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<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)

3
test/adjoint.cpp
View File

@ -143,6 +143,9 @@ template<typename MatrixType> void adjoint(const MatrixType& m)
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>

16
test/triangular.cpp
View File

@ -129,6 +129,22 @@ template<typename MatrixType> void triangular_square(const MatrixType& m)
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());
}