Lazy structured*diagonal evaluator + SelfAdjointView complex-scalar guard libeigen/eigen!2503 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
diff --git a/Eigen/src/Core/ProductEvaluators.h b/Eigen/src/Core/ProductEvaluators.h index 2e4f450..08afb1f 100644 --- a/Eigen/src/Core/ProductEvaluators.h +++ b/Eigen/src/Core/ProductEvaluators.h
@@ -1273,6 +1273,140 @@ #endif }; +// Lazy evaluators for triangular/selfadjoint × diagonal products. +// +// The triangular-assignment dispatcher in TriangularMatrix.h routes +// `triangularView() ?= structured * diagonal` through call_triangular_assignment_loop, +// which constructs an evaluator over the Product expression. Without these +// specializations, evaluator<Product<TriangularView, Diagonal>> falls back to the +// default product_evaluator that materializes a full PlainObject — defeating the +// triangular-aware dispatcher. The kernel reads only the destination's active +// triangle one coefficient at a time, so a coeff()-only evaluator avoids the +// temporary at no loss in functionality. + +template <int Mode, int ProductOrder, typename MatrixType, typename DiagonalType, typename Derived> +struct triangular_diagonal_product_lazy_evaluator_base : evaluator_base<Derived> { + using Scalar = typename ScalarBinaryOpTraits<typename MatrixType::Scalar, typename DiagonalType::Scalar>::ReturnType; + using MatrixScalar = typename MatrixType::Scalar; + + enum { + CoeffReadCost = int(NumTraits<Scalar>::MulCost) + int(evaluator<MatrixType>::CoeffReadCost) + + int(evaluator<DiagonalType>::CoeffReadCost), + Flags = HereditaryBits & static_cast<unsigned int>(evaluator<MatrixType>::Flags), + Alignment = 0 + }; + + EIGEN_DEVICE_FUNC triangular_diagonal_product_lazy_evaluator_base(const MatrixType& mat, const DiagonalType& diag) + : m_diagImpl(diag), m_matImpl(mat) {} + + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar coeff(Index row, Index col) const { + const bool inActive = ((Mode & Upper) == Upper) ? (row <= col) : (row >= col); + if (!inActive) return Scalar(0); + if (row == col) { + if ((Mode & UnitDiag) == UnitDiag) { + return ProductOrder == OnTheLeft ? Scalar(m_diagImpl.coeff(row) * MatrixScalar(1)) + : Scalar(MatrixScalar(1) * m_diagImpl.coeff(col)); + } + if ((Mode & ZeroDiag) == ZeroDiag) return Scalar(0); + } + return ProductOrder == OnTheLeft ? Scalar(m_diagImpl.coeff(row) * m_matImpl.coeff(row, col)) + : Scalar(m_matImpl.coeff(row, col) * m_diagImpl.coeff(col)); + } + + protected: + evaluator<DiagonalType> m_diagImpl; + evaluator<MatrixType> m_matImpl; +}; + +// Triangular × Diagonal +template <typename Lhs, typename Rhs, int ProductKind, int ProductTag> +struct product_evaluator<Product<Lhs, Rhs, ProductKind>, ProductTag, TriangularShape, DiagonalShape> + : triangular_diagonal_product_lazy_evaluator_base< + Lhs::Mode, OnTheRight, typename Lhs::MatrixType, typename Rhs::DiagonalVectorType, + product_evaluator<Product<Lhs, Rhs, ProductKind>, ProductTag, TriangularShape, DiagonalShape>> { + using XprType = Product<Lhs, Rhs, ProductKind>; + using Base = triangular_diagonal_product_lazy_evaluator_base< + Lhs::Mode, OnTheRight, typename Lhs::MatrixType, typename Rhs::DiagonalVectorType, + product_evaluator<XprType, ProductTag, TriangularShape, DiagonalShape>>; + + EIGEN_DEVICE_FUNC explicit product_evaluator(const XprType& xpr) + : Base(xpr.lhs().nestedExpression(), xpr.rhs().diagonal()) {} +}; + +// Diagonal × Triangular +template <typename Lhs, typename Rhs, int ProductKind, int ProductTag> +struct product_evaluator<Product<Lhs, Rhs, ProductKind>, ProductTag, DiagonalShape, TriangularShape> + : triangular_diagonal_product_lazy_evaluator_base< + Rhs::Mode, OnTheLeft, typename Rhs::MatrixType, typename Lhs::DiagonalVectorType, + product_evaluator<Product<Lhs, Rhs, ProductKind>, ProductTag, DiagonalShape, TriangularShape>> { + using XprType = Product<Lhs, Rhs, ProductKind>; + using Base = triangular_diagonal_product_lazy_evaluator_base< + Rhs::Mode, OnTheLeft, typename Rhs::MatrixType, typename Lhs::DiagonalVectorType, + product_evaluator<XprType, ProductTag, DiagonalShape, TriangularShape>>; + + EIGEN_DEVICE_FUNC explicit product_evaluator(const XprType& xpr) + : Base(xpr.rhs().nestedExpression(), xpr.lhs().diagonal()) {} +}; + +// Dense SelfAdjointView statically rejects the Upper|Lower mode (only one half is stored), so the +// off-stored coefficient is always reconstructed by conjugating its mirror. +template <int Mode, int ProductOrder, typename MatrixType, typename DiagonalType, typename Derived> +struct selfadjoint_diagonal_product_lazy_evaluator_base : evaluator_base<Derived> { + using Scalar = typename ScalarBinaryOpTraits<typename MatrixType::Scalar, typename DiagonalType::Scalar>::ReturnType; + + enum { + CoeffReadCost = int(NumTraits<Scalar>::MulCost) + int(evaluator<MatrixType>::CoeffReadCost) + + int(evaluator<DiagonalType>::CoeffReadCost), + Flags = HereditaryBits & static_cast<unsigned int>(evaluator<MatrixType>::Flags), + Alignment = 0 + }; + + EIGEN_DEVICE_FUNC selfadjoint_diagonal_product_lazy_evaluator_base(const MatrixType& mat, const DiagonalType& diag) + : m_diagImpl(diag), m_matImpl(mat) {} + + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar coeff(Index row, Index col) const { + const bool storedHere = ((Mode & Upper) == Upper) ? (row <= col) : (row >= col); + const Scalar matCoeff = + storedHere ? Scalar(m_matImpl.coeff(row, col)) : Scalar(numext::conj(m_matImpl.coeff(col, row))); + return ProductOrder == OnTheLeft ? Scalar(m_diagImpl.coeff(row) * matCoeff) + : Scalar(matCoeff * m_diagImpl.coeff(col)); + } + + protected: + evaluator<DiagonalType> m_diagImpl; + evaluator<MatrixType> m_matImpl; +}; + +// SelfAdjoint × Diagonal +template <typename Lhs, typename Rhs, int ProductKind, int ProductTag> +struct product_evaluator<Product<Lhs, Rhs, ProductKind>, ProductTag, SelfAdjointShape, DiagonalShape> + : selfadjoint_diagonal_product_lazy_evaluator_base< + Lhs::Mode, OnTheRight, typename Lhs::MatrixType, typename Rhs::DiagonalVectorType, + product_evaluator<Product<Lhs, Rhs, ProductKind>, ProductTag, SelfAdjointShape, DiagonalShape>> { + using XprType = Product<Lhs, Rhs, ProductKind>; + using Base = selfadjoint_diagonal_product_lazy_evaluator_base< + Lhs::Mode, OnTheRight, typename Lhs::MatrixType, typename Rhs::DiagonalVectorType, + product_evaluator<XprType, ProductTag, SelfAdjointShape, DiagonalShape>>; + + EIGEN_DEVICE_FUNC explicit product_evaluator(const XprType& xpr) + : Base(xpr.lhs().nestedExpression(), xpr.rhs().diagonal()) {} +}; + +// Diagonal × SelfAdjoint +template <typename Lhs, typename Rhs, int ProductKind, int ProductTag> +struct product_evaluator<Product<Lhs, Rhs, ProductKind>, ProductTag, DiagonalShape, SelfAdjointShape> + : selfadjoint_diagonal_product_lazy_evaluator_base< + Rhs::Mode, OnTheLeft, typename Rhs::MatrixType, typename Lhs::DiagonalVectorType, + product_evaluator<Product<Lhs, Rhs, ProductKind>, ProductTag, DiagonalShape, SelfAdjointShape>> { + using XprType = Product<Lhs, Rhs, ProductKind>; + using Base = selfadjoint_diagonal_product_lazy_evaluator_base< + Rhs::Mode, OnTheLeft, typename Rhs::MatrixType, typename Lhs::DiagonalVectorType, + product_evaluator<XprType, ProductTag, DiagonalShape, SelfAdjointShape>>; + + EIGEN_DEVICE_FUNC explicit product_evaluator(const XprType& xpr) + : Base(xpr.rhs().nestedExpression(), xpr.lhs().diagonal()) {} +}; + /*************************************************************************** * Products with permutation matrices ***************************************************************************/
diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h index 4056fdc..2546d4a 100644 --- a/Eigen/src/Core/SelfAdjointView.h +++ b/Eigen/src/Core/SelfAdjointView.h
@@ -108,12 +108,20 @@ /** \sa MatrixBase::operator*=() */ EIGEN_DEVICE_FUNC SelfAdjointView& operator*=(const Scalar& other) { + eigen_assert(numext::imag(other) == typename NumTraits<Scalar>::Real(0) && + "SelfAdjointView in-place scaling requires a real scalar; " + "scaling only the stored triangle by a non-real scalar would " + "leave conj(other) on the unstored half."); m_matrix.template triangularView<UpLo>() *= other; return *this; } /** \sa DenseBase::operator/=() */ EIGEN_DEVICE_FUNC SelfAdjointView& operator/=(const Scalar& other) { + eigen_assert(numext::imag(other) == typename NumTraits<Scalar>::Real(0) && + "SelfAdjointView in-place division requires a real scalar; " + "dividing only the stored triangle by a non-real scalar would " + "leave conj(other) on the unstored half."); m_matrix.template triangularView<UpLo>() /= other; return *this; }
diff --git a/test/diagonalmatrices.cpp b/test/diagonalmatrices.cpp index e61b4b6..97046de 100644 --- a/test/diagonalmatrices.cpp +++ b/test/diagonalmatrices.cpp
@@ -230,6 +230,28 @@ no_malloc_result.noalias() = dynamic_d.asDiagonal() * dynamic_m.template triangularView<UnitUpper>(); internal::set_is_malloc_allowed(true); VERIFY_IS_APPROX(no_malloc_result, dynamic_expected); + + // Triangular destination assignment from a structured*diagonal product must go through the + // lazy product_evaluator and avoid materializing a full PlainObject temporary. + no_malloc_result = dynamic_m; + dynamic_expected = dynamic_m; + dynamic_expected.template triangularView<Upper>() = + MatrixXd(dynamic_m.template triangularView<Upper>()) * dynamic_d.asDiagonal(); + internal::set_is_malloc_allowed(false); + no_malloc_result.template triangularView<Upper>() = + dynamic_m.template triangularView<Upper>() * dynamic_d.asDiagonal(); + internal::set_is_malloc_allowed(true); + VERIFY_IS_APPROX(no_malloc_result, dynamic_expected); + + no_malloc_result = dynamic_m; + dynamic_expected = dynamic_m; + dynamic_expected.template triangularView<Lower>() = + dynamic_d.asDiagonal() * MatrixXd(dynamic_m.template triangularView<Lower>()); + internal::set_is_malloc_allowed(false); + no_malloc_result.template triangularView<Lower>() = + dynamic_d.asDiagonal() * dynamic_m.template triangularView<Lower>(); + internal::set_is_malloc_allowed(true); + VERIFY_IS_APPROX(no_malloc_result, dynamic_expected); } template <int> @@ -290,6 +312,28 @@ no_malloc_result.noalias() = dynamic_d.asDiagonal() * dynamic_m.template selfadjointView<Upper>(); internal::set_is_malloc_allowed(true); VERIFY_IS_APPROX(no_malloc_result, dynamic_expected); + + // Triangular destination assignment from selfadjoint*diagonal must use the lazy + // product_evaluator (with conjugate-mirror coeffs) and not materialize a temporary. + no_malloc_result = dynamic_m; + dynamic_expected = dynamic_m; + dynamic_expected.template triangularView<Lower>() = + MatrixXcd(dynamic_m.template selfadjointView<Upper>()) * dynamic_d.asDiagonal(); + internal::set_is_malloc_allowed(false); + no_malloc_result.template triangularView<Lower>() = + dynamic_m.template selfadjointView<Upper>() * dynamic_d.asDiagonal(); + internal::set_is_malloc_allowed(true); + VERIFY_IS_APPROX(no_malloc_result, dynamic_expected); + + no_malloc_result = dynamic_m; + dynamic_expected = dynamic_m; + dynamic_expected.template triangularView<Upper>() = + dynamic_d.asDiagonal() * MatrixXcd(dynamic_m.template selfadjointView<Lower>()); + internal::set_is_malloc_allowed(false); + no_malloc_result.template triangularView<Upper>() = + dynamic_d.asDiagonal() * dynamic_m.template selfadjointView<Lower>(); + internal::set_is_malloc_allowed(true); + VERIFY_IS_APPROX(no_malloc_result, dynamic_expected); } EIGEN_DECLARE_TEST(diagonalmatrices) {
diff --git a/test/product_selfadjoint.cpp b/test/product_selfadjoint.cpp index 5eda21d..42cc0d6 100644 --- a/test/product_selfadjoint.cpp +++ b/test/product_selfadjoint.cpp
@@ -12,6 +12,7 @@ template <typename MatrixType> void product_selfadjoint(const MatrixType& m) { typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; typedef Matrix<Scalar, 1, MatrixType::RowsAtCompileTime> RowVectorType; @@ -84,17 +85,20 @@ m3.template triangularView<Lower>() -= m1; VERIFY_IS_APPROX(m2, m3); + // SelfAdjointView in-place scaling/division accept real scalars only: scaling only the + // stored triangle by a non-real factor would leave conj(other) on the unstored half. + const RealScalar real_s1 = numext::real(s1); + const RealScalar real_divisor = numext::real(s2) + RealScalar(2); m2.setRandom(); m3 = m2; - m2.template selfadjointView<Upper>() *= s1; - m3.template triangularView<Upper>() *= s1; + m2.template selfadjointView<Upper>() *= real_s1; + m3.template triangularView<Upper>() *= real_s1; VERIFY_IS_APPROX(m2, m3); m2.setRandom(); m3 = m2; - const Scalar divisor = s2 + Scalar(2); - m2.template selfadjointView<Lower>() /= divisor; - m3.template triangularView<Lower>() /= divisor; + m2.template selfadjointView<Lower>() /= real_divisor; + m3.template triangularView<Lower>() /= real_divisor; VERIFY_IS_APPROX(m2, m3); m2.setRandom();
diff --git a/test/sparse_product.cpp b/test/sparse_product.cpp index f001a90..8ec80e9 100644 --- a/test/sparse_product.cpp +++ b/test/sparse_product.cpp
@@ -370,6 +370,10 @@ VERIFY_IS_APPROX(x = scale.asDiagonal() * mLo.template selfadjointView<Lower>(), refX = scale.asDiagonal() * refS); VERIFY_IS_APPROX(x = mUp.template selfadjointView<Upper>() * scale.asDiagonal(), refX = refS * scale.asDiagonal()); VERIFY_IS_APPROX(x = scale.asDiagonal() * mUp.template selfadjointView<Upper>(), refX = scale.asDiagonal() * refS); + VERIFY_IS_APPROX(x = mS.template selfadjointView<Upper | Lower>() * scale.asDiagonal(), + refX = refS * scale.asDiagonal()); + VERIFY_IS_APPROX(x = scale.asDiagonal() * mS.template selfadjointView<Upper | Lower>(), + refX = scale.asDiagonal() * refS); // sparse selfadjointView with sparse matrices SparseMatrixType mSres(rows, rows); @@ -385,6 +389,10 @@ refX = refS * scale.asDiagonal()); VERIFY_IS_APPROX(mSres = scale.asDiagonal() * mUp.template selfadjointView<Upper>(), refX = scale.asDiagonal() * refS); + VERIFY_IS_APPROX(mSres = mS.template selfadjointView<Upper | Lower>() * scale.asDiagonal(), + refX = refS * scale.asDiagonal()); + VERIFY_IS_APPROX(mSres = scale.asDiagonal() * mS.template selfadjointView<Upper | Lower>(), + refX = scale.asDiagonal() * refS); // sparse triangularView with dense matrices VERIFY_IS_APPROX(x = mA.template triangularView<Upper>() * b, refX = refA.template triangularView<Upper>() * b);