diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h index 62d072946..ff450eb28 100644 --- a/Eigen/src/Core/SelfAdjointView.h +++ b/Eigen/src/Core/SelfAdjointView.h @@ -114,6 +114,12 @@ class SelfAdjointView : public TriangularBase return Product(lhs.derived(), rhs); } + EIGEN_DEVICE_FUNC const + SelfAdjointView + operator*(const Scalar& s) const { + return (nestedExpression() * s).template selfadjointView(); + } + friend EIGEN_DEVICE_FUNC const SelfAdjointView operator*(const Scalar& s, const SelfAdjointView& mat) { diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index e219e51e3..2b0f56bea 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -413,6 +413,21 @@ class TriangularViewImpl : public TriangularBase(lhs.derived(), rhs.derived()); } + // Scaling a unit triangular view would break its implicit unit diagonal, so only non-unit modes participate. + template = 0> + EIGEN_DEVICE_FUNC const + TriangularView + operator*(const Scalar& s) const { + return (derived().nestedExpression() * s).template triangularView(); + } + + template = 0> + friend EIGEN_DEVICE_FUNC const + TriangularView + operator*(const Scalar& s, const TriangularViewImpl& mat) { + return (s * mat.derived().nestedExpression()).template triangularView(); + } + /** \returns the product of the inverse of \c *this with \a other, \a *this being triangular. * * This function computes the inverse-matrix matrix product inverse(\c *this) * \a other if diff --git a/Eigen/src/SparseCore/SparseSelfAdjointView.h b/Eigen/src/SparseCore/SparseSelfAdjointView.h index 9b290ddd8..32fbcd4fb 100644 --- a/Eigen/src/SparseCore/SparseSelfAdjointView.h +++ b/Eigen/src/SparseCore/SparseSelfAdjointView.h @@ -64,6 +64,8 @@ class SparseSelfAdjointView : public EigenBase VectorI; typedef typename internal::ref_selector::non_const_type MatrixTypeNested; typedef internal::remove_all_t MatrixTypeNested_; + typedef SparseMatrix + PlainObject; explicit inline SparseSelfAdjointView(MatrixType& matrix) : m_matrix(matrix) { eigen_assert(rows() == cols() && "SelfAdjointView is only for squared matrices"); @@ -114,6 +116,16 @@ class SparseSelfAdjointView : public EigenBase(lhs.derived(), rhs); } + // Scalar multiplication intentionally materializes the full matrix, unlike dense SelfAdjointView's lazy wrapper, + // matching the existing SparseSelfAdjointView products. + PlainObject operator*(const Scalar& s) const { return s * *this; } + + friend PlainObject operator*(const Scalar& s, const SparseSelfAdjointView& mat) { + PlainObject res(mat); + res *= s; + return res; + } + /** Perform a symmetric rank K update of the selfadjoint matrix \c *this: * \f$ this = this + \alpha ( u u^* ) \f$ where \a u is a vector or matrix. * diff --git a/test/selfadjoint.cpp b/test/selfadjoint.cpp index 01e3806dc..65d4e7982 100644 --- a/test/selfadjoint.cpp +++ b/test/selfadjoint.cpp @@ -42,6 +42,18 @@ void selfadjoint(const MatrixType& m) { m4 = m2; m4 -= m1.template selfadjointView(); VERIFY_IS_APPROX(m4, m2 - m3); + + Scalar s = internal::random(); + + m4 = s * m1.template selfadjointView(); + VERIFY_IS_APPROX(m4, MatrixType((s * m1).template selfadjointView())); + m4 = m1.template selfadjointView() * s; + VERIFY_IS_APPROX(m4, MatrixType((m1 * s).template selfadjointView())); + + m4 = s * m1.template selfadjointView(); + VERIFY_IS_APPROX(m4, MatrixType((s * m1).template selfadjointView())); + m4 = m1.template selfadjointView() * s; + VERIFY_IS_APPROX(m4, MatrixType((m1 * s).template selfadjointView())); } void bug_159() { diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp index 4f1d4473e..d583d9d13 100644 --- a/test/sparse_basic.cpp +++ b/test/sparse_basic.cpp @@ -796,10 +796,28 @@ void sparse_basic(const SparseMatrixType& ref) { m3 -= m2.template selfadjointView(); VERIFY_IS_APPROX(m3, refMat3); + Scalar s2 = internal::random(); + refMat3 = DenseMatrix(refMat2.template selfadjointView()); + refMat3 *= s2; + SparseMatrixType m4 = s2 * m2.template selfadjointView(); + VERIFY_IS_APPROX(m4, refMat3); + refMat3 = DenseMatrix(refMat2.template selfadjointView()); + refMat3 *= s2; + m4 = m2.template selfadjointView() * s2; + VERIFY_IS_APPROX(m4, refMat3); + refMat3 = DenseMatrix(refMat2.template selfadjointView()); + refMat3 *= s2; + m4 = s2 * m2.template selfadjointView(); + VERIFY_IS_APPROX(m4, refMat3); + refMat3 = DenseMatrix(refMat2.template selfadjointView()); + refMat3 *= s2; + m4 = m2.template selfadjointView() * s2; + VERIFY_IS_APPROX(m4, refMat3); + // selfadjointView only works for square matrices: - SparseMatrixType m4(rows, rows + 1); - VERIFY_RAISES_ASSERT(m4.template selfadjointView()); - VERIFY_RAISES_ASSERT(m4.template selfadjointView()); + SparseMatrixType m5(rows, rows + 1); + VERIFY_RAISES_ASSERT(m5.template selfadjointView()); + VERIFY_RAISES_ASSERT(m5.template selfadjointView()); } // test sparseView diff --git a/test/triangular.cpp b/test/triangular.cpp index a539715f5..ea4420b42 100644 --- a/test/triangular.cpp +++ b/test/triangular.cpp @@ -13,6 +13,55 @@ #include "main.h" +template +struct has_left_scalar_multiply : std::false_type {}; + +template +struct has_left_scalar_multiply< + ViewType, internal::void_t() * std::declval())>> + : std::true_type {}; + +template +struct has_right_scalar_multiply : std::false_type {}; + +template +struct has_right_scalar_multiply< + ViewType, internal::void_t() * std::declval())>> + : std::true_type {}; + +template +void triangular_scalar_multiply(const MatrixType& m) { + typedef typename MatrixType::Scalar Scalar; + + const Index rows = m.rows(); + const Index cols = m.cols(); + + const Scalar s = internal::random(); + const MatrixType triangular = MatrixType::Random(rows, cols); + + VERIFY_IS_APPROX((s * triangular.template triangularView()).toDenseMatrix(), + (s * triangular).template triangularView().toDenseMatrix()); + VERIFY_IS_APPROX((triangular.template triangularView() * s).toDenseMatrix(), + (triangular * s).template triangularView().toDenseMatrix()); +} + +template +void triangular_scalar_multiply_sfinae() { + typedef decltype(std::declval().template triangularView()) LowerView; + typedef decltype(std::declval().template triangularView()) StrictlyLowerView; + typedef decltype(std::declval().template triangularView()) StrictlyUpperView; + typedef decltype(std::declval().template triangularView()) UnitLowerView; + + STATIC_CHECK((has_left_scalar_multiply::value)); + STATIC_CHECK((has_right_scalar_multiply::value)); + STATIC_CHECK((has_left_scalar_multiply::value)); + STATIC_CHECK((has_right_scalar_multiply::value)); + STATIC_CHECK((has_left_scalar_multiply::value)); + STATIC_CHECK((has_right_scalar_multiply::value)); + STATIC_CHECK((!has_left_scalar_multiply::value)); + STATIC_CHECK((!has_right_scalar_multiply::value)); +} + template void triangular_deprecated(const MatrixType& m) { Index rows = m.rows(); @@ -42,6 +91,8 @@ void triangular_square(const MatrixType& m) { typedef typename NumTraits::Real RealScalar; typedef Matrix VectorType; + triangular_scalar_multiply_sfinae(); + RealScalar largerEps = 10 * test_precision(); Index rows = m.rows(); @@ -151,6 +202,10 @@ void triangular_square(const MatrixType& m) { m6.setRandom(); VERIFY_IS_APPROX(m1.template triangularView() * m5, m3 * m5); VERIFY_IS_APPROX(m6 * m1.template triangularView(), m6 * m3); + triangular_scalar_multiply(m1); + triangular_scalar_multiply(m1); + triangular_scalar_multiply(m1); + triangular_scalar_multiply(m1); m1up = m1.template triangularView(); VERIFY_IS_APPROX(m1.template selfadjointView().template triangularView().toDenseMatrix(), m1up); @@ -228,6 +283,10 @@ void triangular_rect(const MatrixType& m) { m1.setZero(); m1.template triangularView() = 3 * m2; VERIFY_IS_APPROX(m3.template triangularView().toDenseMatrix(), m1); + triangular_scalar_multiply(m1); + triangular_scalar_multiply(m1); + triangular_scalar_multiply(m1); + triangular_scalar_multiply(m1); m1.setRandom(); m2 = m1.template triangularView(); VERIFY(m2.isUpperTriangular());