From 0cb4f32e12a190f43f794a8258661bade74f2eb2 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 22 Jul 2009 23:12:22 +0200 Subject: [PATCH] implement high level API for SYMM and fix a couple of bugs related to complex --- Eigen/src/Core/MatrixBase.h | 4 + Eigen/src/Core/ReturnByValue.h | 36 +++++++- Eigen/src/Core/SelfAdjointView.h | 89 +++++++++++++++++-- Eigen/src/Core/products/GeneralMatrixMatrix.h | 7 +- .../Core/products/SelfadjointMatrixMatrix.h | 49 +++++----- Eigen/src/Core/util/BlasUtil.h | 2 +- Eigen/src/Core/util/Macros.h | 1 + test/product_selfadjoint.cpp | 68 ++++++++++---- 8 files changed, 200 insertions(+), 56 deletions(-) diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 0a3342758..3df518d1a 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -300,6 +300,10 @@ template class MatrixBase template Derived& operator=(const ReturnByValue& func); + template + Derived& operator+=(const ReturnByValue& func); + template + Derived& operator-=(const ReturnByValue& func); #ifndef EIGEN_PARSED_BY_DOXYGEN /** Copies \a other into *this without evaluating other. \returns a reference to *this. */ diff --git a/Eigen/src/Core/ReturnByValue.h b/Eigen/src/Core/ReturnByValue.h index fc9211cdd..72f4a8f43 100644 --- a/Eigen/src/Core/ReturnByValue.h +++ b/Eigen/src/Core/ReturnByValue.h @@ -46,20 +46,36 @@ struct ei_nested >, n, EvalTyp template class ReturnByValue { public: - template - inline void evalTo(Dest& dst) const + template inline void evalTo(Dest& dst) const { static_cast(this)->evalTo(dst); } + template inline void addTo(Dest& dst) const + { static_cast(this)->_addTo(dst); } + template inline void subTo(Dest& dst) const + { static_cast(this)->_subTo(dst); } + template inline void _addTo(Dest& dst) const + { EvalType res; evalTo(res); dst += res; } + template inline void _subTo(Dest& dst) const + { EvalType res; evalTo(res); dst -= res; } }; template class ReturnByValue > : public MatrixBase > > { + typedef Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> EvalType; public: EIGEN_GENERIC_PUBLIC_INTERFACE(ReturnByValue) template inline void evalTo(Dest& dst) const { static_cast(this)->evalTo(dst); } + template inline void addTo(Dest& dst) const + { static_cast(this)->_addTo(dst); } + template inline void subTo(Dest& dst) const + { static_cast(this)->_subTo(dst); } + template inline void _addTo(Dest& dst) const + { EvalType res; evalTo(res); dst += res; } + template inline void _subTo(Dest& dst) const + { EvalType res; evalTo(res); dst -= res; } }; template @@ -70,4 +86,20 @@ Derived& MatrixBase::operator=(const ReturnByValue +template +Derived& MatrixBase::operator+=(const ReturnByValue& other) +{ + other.addTo(derived()); + return derived(); +} + +template +template +Derived& MatrixBase::operator-=(const ReturnByValue& other) +{ + other.subTo(derived()); + return derived(); +} + #endif // EIGEN_RETURNBYVALUE_H diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h index a8ea7cd62..7f5fd7533 100644 --- a/Eigen/src/Core/SelfAdjointView.h +++ b/Eigen/src/Core/SelfAdjointView.h @@ -53,8 +53,8 @@ struct ei_traits > : ei_traits -struct ei_selfadjoint_vector_product_returntype; +template +struct ei_selfadjoint_matrix_product_returntype; // FIXME could also be called SelfAdjointWrapper to be consistent with DiagonalWrapper ?? template class SelfAdjointView @@ -97,13 +97,12 @@ template class SelfAdjointView /** \internal */ const MatrixType& _expression() const { return m_matrix; } - /** Efficient self-adjoint matrix times vector product */ - // TODO this product is far to be ready + /** Efficient self-adjoint matrix times vector/matrix product */ template - ei_selfadjoint_vector_product_returntype + ei_selfadjoint_matrix_product_returntype operator*(const MatrixBase& rhs) const { - return ei_selfadjoint_vector_product_returntype(*this, rhs.derived()); + return ei_selfadjoint_matrix_product_returntype(*this, rhs.derived()); } /** Perform a symmetric rank 2 update of the selfadjoint matrix \c *this: @@ -165,13 +164,13 @@ struct ei_triangular_assignment_selector -struct ei_selfadjoint_vector_product_returntype - : public ReturnByValue, +struct ei_selfadjoint_matrix_product_returntype + : public ReturnByValue, Matrix::Scalar, Rhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> > { typedef typename ei_cleantype::type RhsNested; - ei_selfadjoint_vector_product_returntype(const Lhs& lhs, const Rhs& rhs) + ei_selfadjoint_matrix_product_returntype(const Lhs& lhs, const Rhs& rhs) : m_lhs(lhs), m_rhs(rhs) {} @@ -194,6 +193,78 @@ struct ei_selfadjoint_vector_product_returntype const typename Rhs::Nested m_rhs; }; +/*************************************************************************** +* Wrapper to ei_product_selfadjoint_matrix +***************************************************************************/ + +template +struct ei_selfadjoint_matrix_product_returntype + : public ReturnByValue, + Matrix::Scalar, + Rhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> > +{ + ei_selfadjoint_matrix_product_returntype(const Lhs& lhs, const Rhs& rhs) + : m_lhs(lhs), m_rhs(rhs) + {} + + typedef typename Lhs::Scalar Scalar; + typedef typename Rhs::Nested RhsNested; + typedef typename ei_cleantype::type _RhsNested; + + typedef typename ei_traits::ExpressionType LhsExpr; + typedef typename LhsExpr::Nested LhsNested; + typedef typename ei_cleantype::type _LhsNested; + + enum { UpLo = ei_traits::Mode&(UpperTriangularBit|LowerTriangularBit) }; + + template inline void _addTo(Dest& dst) const + { evalTo(dst,1); } + template inline void _subTo(Dest& dst) const + { evalTo(dst,-1); } + + template void evalTo(Dest& dst) const + { + dst.resize(m_lhs.rows(), m_rhs.cols()); + dst.setZero(); + evalTo(dst,1); + } + + template 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::type _ActualLhsType; + typedef typename ei_cleantype::type _ActualRhsType; + + const ActualLhsType lhs = LhsBlasTraits::extract(m_lhs._expression()); + const ActualRhsType rhs = RhsBlasTraits::extract(m_rhs); + + Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(m_lhs._expression()) + * RhsBlasTraits::extractScalarFactor(m_rhs); + + ei_product_selfadjoint_matrix::Flags &RowMajorBit) ? RowMajor : ColMajor, true, + NumTraits::IsComplex && EIGEN_LOGICAL_XOR(UpLo==UpperTriangular,bool(LhsBlasTraits::NeedToConjugate)), + ei_traits::Flags &RowMajorBit ? RowMajor : ColMajor, false, bool(RhsBlasTraits::NeedToConjugate), + ei_traits::Flags&RowMajorBit ? RowMajor : ColMajor> + ::run( + lhs.rows(), rhs.cols(), // sizes + &lhs.coeff(0,0), lhs.stride(), // lhs info + &rhs.coeff(0,0), rhs.stride(), // rhs info + &dst.coeffRef(0,0), dst.stride(), // result info + actualAlpha // alpha + ); + } + + const Lhs m_lhs; + const RhsNested m_rhs; +}; + /*************************************************************************** * Implementation of MatrixBase methods ***************************************************************************/ diff --git a/Eigen/src/Core/products/GeneralMatrixMatrix.h b/Eigen/src/Core/products/GeneralMatrixMatrix.h index 1c48a5ed4..b2f51ca5b 100644 --- a/Eigen/src/Core/products/GeneralMatrixMatrix.h +++ b/Eigen/src/Core/products/GeneralMatrixMatrix.h @@ -334,22 +334,23 @@ struct ei_gebp_kernel }; // pack a block of the lhs -template +template struct ei_gemm_pack_lhs { void operator()(Scalar* blockA, const EIGEN_RESTRICT Scalar* _lhs, int lhsStride, int actual_kc, int actual_mc) { + ei_conj_if::IsComplex && Conjugate> cj; ei_const_blas_data_mapper lhs(_lhs,lhsStride); int count = 0; const int peeled_mc = (actual_mc/mr)*mr; for(int i=0; i rhs(_rhs,rhsStride); - // first part: standard case + // first part: normal case for(int j2=0; j2 - ::run(rows, cols, rhs, rhsStride, lhs, lhsStride, res, resStride, alpha); + EIGEN_LOGICAL_XOR(LhsSelfAdjoint,LhsStorageOrder==RowMajor) ? ColMajor : RowMajor, + LhsSelfAdjoint, NumTraits::IsComplex && !ConjugateLhs, ColMajor> + ::run(cols, rows, rhs, rhsStride, lhs, lhsStride, res, resStride, alpha); } }; @@ -254,8 +253,8 @@ struct ei_product_selfadjoint_matrix() + // transposed packed copy + ei_gemm_pack_lhs() (blockA, &lhs(k2, i2), lhsStride, actual_kc, actual_mc); gebp_kernel(res, resStride, blockA, blockB, actual_mc, actual_kc, packet_cols, i2, cols); @@ -273,7 +272,7 @@ struct ei_product_selfadjoint_matrix() + ei_gemm_pack_lhs() (blockA, &lhs(i2, k2), lhsStride, actual_kc, actual_mc); gebp_kernel(res, resStride, blockA, blockB, actual_mc, actual_kc, packet_cols, i2, cols); diff --git a/Eigen/src/Core/util/BlasUtil.h b/Eigen/src/Core/util/BlasUtil.h index 25829652f..10cbfb069 100644 --- a/Eigen/src/Core/util/BlasUtil.h +++ b/Eigen/src/Core/util/BlasUtil.h @@ -35,7 +35,7 @@ struct ei_gebp_kernel; template struct ei_gemm_pack_rhs; -template +template struct ei_gemm_pack_lhs; template< diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h index c9f2f4d40..956b609a0 100644 --- a/Eigen/src/Core/util/Macros.h +++ b/Eigen/src/Core/util/Macros.h @@ -277,5 +277,6 @@ _EIGEN_GENERIC_PUBLIC_INTERFACE(Derived, Eigen::MatrixBase) #define EIGEN_ENUM_MIN(a,b) (((int)a <= (int)b) ? (int)a : (int)b) #define EIGEN_ENUM_MAX(a,b) (((int)a >= (int)b) ? (int)a : (int)b) +#define EIGEN_LOGICAL_XOR(a,b) (((a) || (b)) && !((a) && (b))) #endif // EIGEN_MACROS_H diff --git a/test/product_selfadjoint.cpp b/test/product_selfadjoint.cpp index 814672696..eca1daacf 100644 --- a/test/product_selfadjoint.cpp +++ b/test/product_selfadjoint.cpp @@ -87,6 +87,53 @@ template void product_selfadjoint(const MatrixType& m) } } +template void symm(const MatrixType& m) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef Matrix Rhs1; + typedef Matrix Rhs2; + typedef Matrix Rhs3; + + int rows = m.rows(); + int cols = m.cols(); + + MatrixType m1 = MatrixType::Random(rows, cols), + m2 = MatrixType::Random(rows, cols); + + m1 = (m1+m1.adjoint()).eval(); + + Rhs1 rhs1 = Rhs1::Random(cols, ei_random(1,320)), rhs12, rhs13; + Rhs2 rhs2 = Rhs2::Random(ei_random(1,320), rows), rhs22, rhs23; + Rhs3 rhs3 = Rhs3::Random(cols, ei_random(1,320)), rhs32, rhs33; + +// Scalar s1 = ei_random(), +// s2 = ei_random(); + + m2 = m1.template triangularView(); + VERIFY_IS_APPROX(rhs12 = m2.template selfadjointView() * rhs1, rhs13 = m1 * rhs1); + + m2 = m1.template triangularView(); + VERIFY_IS_APPROX(rhs12 = m2.template selfadjointView() * rhs1, rhs13 = m1 * rhs1); + + m2 = m1.template triangularView(); + VERIFY_IS_APPROX(rhs22 = m2.template selfadjointView() * rhs2.adjoint(), rhs23 = m1 * rhs2.adjoint()); + + m2 = m1.template triangularView(); + VERIFY_IS_APPROX(rhs22 = m2.template selfadjointView() * rhs2.adjoint(), rhs23 = m1 * rhs2.adjoint()); + + m2 = m1.template triangularView(); + VERIFY_IS_APPROX(rhs22 = m2.adjoint().template selfadjointView() * rhs2.adjoint(), + rhs23 = m1.adjoint() * rhs2.adjoint()); + + // test row major = <...> + m2 = m1.template triangularView(); + VERIFY_IS_APPROX(rhs32 = m2.template selfadjointView() * rhs3, rhs33 = m1 * rhs3); + + m2 = m1.template triangularView(); + VERIFY_IS_APPROX(rhs32 = m2.adjoint().template selfadjointView() * rhs3.conjugate(), + rhs33 = m1.adjoint() * rhs3.conjugate()); +} void test_product_selfadjoint() { for(int i = 0; i < g_repeat ; i++) { @@ -102,21 +149,10 @@ void test_product_selfadjoint() for(int i = 0; i < g_repeat ; i++) { - int size = ei_random(10,1024); - int cols = ei_random(10,320); - MatrixXf A = MatrixXf::Random(size,size); - MatrixXf B = MatrixXf::Random(size,cols); - MatrixXf C = MatrixXf::Random(size,cols); - MatrixXf R = MatrixXf::Random(size,cols); - A = (A+A.transpose()).eval(); - - R = C + (A * B).eval(); - - A.corner(TopRight,size-1,size-1).triangularView().setZero(); - - ei_product_selfadjoint_matrix - (size, A.data(), A.stride(), B.data(), B.stride(), false, B.cols(), C.data(), C.stride(), 1); -// std::cerr << A << "\n\n" << C << "\n\n" << R << "\n\n"; - VERIFY_IS_APPROX(C,R); + int s; + s = ei_random(10,320); + CALL_SUBTEST( symm(MatrixXf(s, s)) ); + s = ei_random(10,320); + CALL_SUBTEST( symm(MatrixXcd(s, s)) ); } }