Unify triangular and self-adjoint view interfaces libeigen/eigen!2486 Closes #784 and #2013
diff --git a/Eigen/src/Core/ProductEvaluators.h b/Eigen/src/Core/ProductEvaluators.h index d789f75..2e4f450 100644 --- a/Eigen/src/Core/ProductEvaluators.h +++ b/Eigen/src/Core/ProductEvaluators.h
@@ -847,6 +847,76 @@ template <int Mode, bool LhsIsTriangular, typename Lhs, bool LhsIsVector, typename Rhs, bool RhsIsVector> struct triangular_product_impl; +template <int ProductOrder> +struct diagonal_product_segment_impl; + +template <> +struct diagonal_product_segment_impl<OnTheLeft> { + template <typename DstSegment, typename Coeffs, typename DiagonalType, typename Alpha> + static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(DstSegment& dst, const Coeffs& coeffs, + const DiagonalType& diagonal, Index begin, Index /*col*/, + const Alpha& alpha) { + dst += alpha * diagonal.segment(begin, dst.size()).cwiseProduct(coeffs); + } +}; + +template <> +struct diagonal_product_segment_impl<OnTheRight> { + template <typename DstSegment, typename Coeffs, typename DiagonalType, typename Alpha> + static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(DstSegment& dst, const Coeffs& coeffs, + const DiagonalType& diagonal, Index /*begin*/, Index col, + const Alpha& alpha) { + dst += alpha * (coeffs * diagonal.coeff(col)); + } +}; + +template <int Mode, int ProductOrder, typename MatrixType, typename DiagonalType> +struct triangular_diagonal_product_impl { + typedef typename MatrixType::Scalar MatrixScalar; + + template <typename Dest, typename Alpha> + static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(Dest& dst, const MatrixType& matrix, + const DiagonalType& diagonal, const Alpha& alpha) { + eigen_assert(((ProductOrder == OnTheLeft && diagonal.size() == matrix.rows()) || + (ProductOrder == OnTheRight && diagonal.size() == matrix.cols())) && + "invalid matrix product"); + + const Index rows = matrix.rows(); + const Index cols = matrix.cols(); + for (Index col = 0; col < cols; ++col) { + if ((Mode & Upper) == Upper) { + const Index end = (std::min)(rows, ((Mode & (UnitDiag | ZeroDiag)) ? col : col + 1)); + addStoredSegment(dst, matrix, diagonal, 0, end, col, alpha); + } else { + const Index begin = ((Mode & (UnitDiag | ZeroDiag)) ? col + 1 : col); + addStoredSegment(dst, matrix, diagonal, begin, rows - begin, col, alpha); + } + + if ((Mode & UnitDiag) == UnitDiag && col < rows) addUnitCoeff(dst, diagonal, col, alpha); + } + } + + private: + template <typename Dest, typename Alpha> + static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void addStoredSegment(Dest& dst, const MatrixType& matrix, + const DiagonalType& diagonal, Index begin, + Index size, Index col, const Alpha& alpha) { + if (size <= 0) return; + auto dstSegment = dst.col(col).segment(begin, size); + diagonal_product_segment_impl<ProductOrder>::run(dstSegment, matrix.col(col).segment(begin, size), diagonal, begin, + col, alpha); + } + + template <typename Dest, typename Alpha> + static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void addUnitCoeff(Dest& dst, const DiagonalType& diagonal, Index index, + const Alpha& alpha) { + if (ProductOrder == OnTheLeft) + dst.coeffRef(index, index) += alpha * (diagonal.coeff(index) * MatrixScalar(1)); + else + dst.coeffRef(index, index) += alpha * (MatrixScalar(1) * diagonal.coeff(index)); + } +}; + template <typename Lhs, typename Rhs, int ProductTag> struct generic_product_impl<Lhs, Rhs, TriangularShape, DenseShape, ProductTag> : generic_product_impl_base<Lhs, Rhs, generic_product_impl<Lhs, Rhs, TriangularShape, DenseShape, ProductTag>> { @@ -871,12 +941,82 @@ } }; +template <typename Lhs, typename Rhs, int ProductTag> +struct generic_product_impl<Lhs, Rhs, TriangularShape, DiagonalShape, ProductTag> + : generic_product_impl_base<Lhs, Rhs, generic_product_impl<Lhs, Rhs, TriangularShape, DiagonalShape, ProductTag>> { + typedef typename Product<Lhs, Rhs>::Scalar Scalar; + + template <typename Dest> + static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, + const Scalar& alpha) { + triangular_diagonal_product_impl<Lhs::Mode, OnTheRight, typename Lhs::MatrixType, + typename Rhs::DiagonalVectorType>::run(dst, lhs.nestedExpression(), rhs.diagonal(), + alpha); + } +}; + +template <typename Lhs, typename Rhs, int ProductTag> +struct generic_product_impl<Lhs, Rhs, DiagonalShape, TriangularShape, ProductTag> + : generic_product_impl_base<Lhs, Rhs, generic_product_impl<Lhs, Rhs, DiagonalShape, TriangularShape, ProductTag>> { + typedef typename Product<Lhs, Rhs>::Scalar Scalar; + + template <typename Dest> + static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, + const Scalar& alpha) { + triangular_diagonal_product_impl<Rhs::Mode, OnTheLeft, typename Rhs::MatrixType, + typename Lhs::DiagonalVectorType>::run(dst, rhs.nestedExpression(), lhs.diagonal(), + alpha); + } +}; + /*************************************************************************** * SelfAdjoint products ***************************************************************************/ template <typename Lhs, int LhsMode, bool LhsIsVector, typename Rhs, int RhsMode, bool RhsIsVector> struct selfadjoint_product_impl; +template <int Mode, int ProductOrder, typename MatrixType, typename DiagonalType> +struct selfadjoint_diagonal_product_impl { + template <typename Dest, typename Alpha> + static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(Dest& dst, const MatrixType& matrix, + const DiagonalType& diagonal, const Alpha& alpha) { + eigen_assert(matrix.rows() == matrix.cols() && "SelfAdjointView is only for squared matrices"); + eigen_assert(diagonal.size() == matrix.rows() && "invalid matrix product"); + + const Index size = matrix.rows(); + for (Index col = 0; col < size; ++col) { + if ((Mode & Upper) == Upper) { + addStoredSegment(dst, matrix, diagonal, 0, col + 1, col, alpha); + addConjugateSegment(dst, matrix, diagonal, col + 1, size - col - 1, col, alpha); + } else { + addConjugateSegment(dst, matrix, diagonal, 0, col, col, alpha); + addStoredSegment(dst, matrix, diagonal, col, size - col, col, alpha); + } + } + } + + private: + template <typename Dest, typename Alpha> + static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void addStoredSegment(Dest& dst, const MatrixType& matrix, + const DiagonalType& diagonal, Index begin, + Index size, Index col, const Alpha& alpha) { + if (size <= 0) return; + auto dstSegment = dst.col(col).segment(begin, size); + diagonal_product_segment_impl<ProductOrder>::run(dstSegment, matrix.col(col).segment(begin, size), diagonal, begin, + col, alpha); + } + + template <typename Dest, typename Alpha> + static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void addConjugateSegment(Dest& dst, const MatrixType& matrix, + const DiagonalType& diagonal, Index begin, + Index size, Index col, const Alpha& alpha) { + if (size <= 0) return; + auto dstSegment = dst.col(col).segment(begin, size); + diagonal_product_segment_impl<ProductOrder>::run( + dstSegment, matrix.row(col).segment(begin, size).conjugate().transpose(), diagonal, begin, col, alpha); + } +}; + template <typename Lhs, typename Rhs, int ProductTag> struct generic_product_impl<Lhs, Rhs, SelfAdjointShape, DenseShape, ProductTag> : generic_product_impl_base<Lhs, Rhs, generic_product_impl<Lhs, Rhs, SelfAdjointShape, DenseShape, ProductTag>> { @@ -901,6 +1041,34 @@ } }; +template <typename Lhs, typename Rhs, int ProductTag> +struct generic_product_impl<Lhs, Rhs, SelfAdjointShape, DiagonalShape, ProductTag> + : generic_product_impl_base<Lhs, Rhs, generic_product_impl<Lhs, Rhs, SelfAdjointShape, DiagonalShape, ProductTag>> { + typedef typename Product<Lhs, Rhs>::Scalar Scalar; + + template <typename Dest> + static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, + const Scalar& alpha) { + selfadjoint_diagonal_product_impl<Lhs::Mode, OnTheRight, typename Lhs::MatrixType, + typename Rhs::DiagonalVectorType>::run(dst, lhs.nestedExpression(), + rhs.diagonal(), alpha); + } +}; + +template <typename Lhs, typename Rhs, int ProductTag> +struct generic_product_impl<Lhs, Rhs, DiagonalShape, SelfAdjointShape, ProductTag> + : generic_product_impl_base<Lhs, Rhs, generic_product_impl<Lhs, Rhs, DiagonalShape, SelfAdjointShape, ProductTag>> { + typedef typename Product<Lhs, Rhs>::Scalar Scalar; + + template <typename Dest> + static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, + const Scalar& alpha) { + selfadjoint_diagonal_product_impl<Rhs::Mode, OnTheLeft, typename Rhs::MatrixType, + typename Lhs::DiagonalVectorType>::run(dst, rhs.nestedExpression(), + lhs.diagonal(), alpha); + } +}; + /*************************************************************************** * Diagonal products ***************************************************************************/
diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h index 59731b8..4056fdc 100644 --- a/Eigen/src/Core/SelfAdjointView.h +++ b/Eigen/src/Core/SelfAdjointView.h
@@ -45,6 +45,7 @@ (~(PacketAccessBit | DirectAccessBit | LinearAccessBit)) // FIXME these flags should be preserved }; }; + } // namespace internal template <typename MatrixType_, unsigned int UpLo> @@ -61,8 +62,6 @@ /** \brief The type of coefficients in this matrix */ typedef typename internal::traits<SelfAdjointView>::Scalar Scalar; typedef typename MatrixType::StorageIndex StorageIndex; - typedef internal::remove_all_t<typename MatrixType::ConjugateReturnType> MatrixConjugateReturnType; - typedef SelfAdjointView<std::add_const_t<MatrixType>, UpLo> ConstSelfAdjointView; enum { Mode = internal::traits<SelfAdjointView>::Mode, @@ -72,47 +71,58 @@ typedef typename MatrixType::PlainObject PlainObject; EIGEN_DEVICE_FUNC explicit inline SelfAdjointView(MatrixType& matrix) : m_matrix(matrix) {} + using Base::operator*; + EIGEN_DEFAULT_COPY_CONSTRUCTOR(SelfAdjointView) - EIGEN_DEVICE_FUNC constexpr Index rows() const noexcept { return m_matrix.rows(); } - EIGEN_DEVICE_FUNC constexpr Index cols() const noexcept { return m_matrix.cols(); } - EIGEN_DEVICE_FUNC constexpr Index outerStride() const noexcept { return m_matrix.outerStride(); } - EIGEN_DEVICE_FUNC constexpr Index innerStride() const noexcept { return m_matrix.innerStride(); } - - /** \sa MatrixBase::coeff() - * \warning the coordinates must fit into the referenced triangular part - */ - EIGEN_DEVICE_FUNC inline Scalar coeff(Index row, Index col) const { - Base::check_coordinates_internal(row, col); - return m_matrix.coeff(row, col); + /** Assigns a matrix expression to the referenced triangular part of the selfadjoint matrix. */ + template <typename OtherDerived> + EIGEN_DEVICE_FUNC SelfAdjointView& operator=(const MatrixBase<OtherDerived>& other) { + m_matrix.template triangularView<UpLo>() = other; + return *this; } - /** \sa MatrixBase::coeffRef() - * \warning the coordinates must fit into the referenced triangular part - */ - EIGEN_DEVICE_FUNC inline Scalar& coeffRef(Index row, Index col) { - EIGEN_STATIC_ASSERT_LVALUE(SelfAdjointView); - Base::check_coordinates_internal(row, col); - return m_matrix.coeffRef(row, col); + /** Assigns a triangular or selfadjoint expression without materializing a dense temporary. */ + template <typename OtherDerived> + EIGEN_DEVICE_FUNC SelfAdjointView& operator=(const TriangularBase<OtherDerived>& other) { + other.evalToLazy(m_matrix); + return *this; + } + + EIGEN_DEVICE_FUNC SelfAdjointView& operator=(const SelfAdjointView& other) { + return *this = static_cast<const Base&>(other); + } + + /** \sa MatrixBase::operator+=() */ + template <typename OtherDerived> + EIGEN_DEVICE_FUNC SelfAdjointView& operator+=(const DenseBase<OtherDerived>& other) { + m_matrix.template triangularView<UpLo>() += other; + return *this; + } + + /** \sa MatrixBase::operator-=() */ + template <typename OtherDerived> + EIGEN_DEVICE_FUNC SelfAdjointView& operator-=(const DenseBase<OtherDerived>& other) { + m_matrix.template triangularView<UpLo>() -= other; + return *this; + } + + /** \sa MatrixBase::operator*=() */ + EIGEN_DEVICE_FUNC SelfAdjointView& operator*=(const Scalar& other) { + m_matrix.template triangularView<UpLo>() *= other; + return *this; + } + + /** \sa DenseBase::operator/=() */ + EIGEN_DEVICE_FUNC SelfAdjointView& operator/=(const Scalar& other) { + m_matrix.template triangularView<UpLo>() /= other; + return *this; } /** \internal */ - EIGEN_DEVICE_FUNC const MatrixTypeNestedCleaned& _expression() const { return m_matrix; } + EIGEN_DEVICE_FUNC constexpr const MatrixTypeNestedCleaned& _expression() const noexcept { return m_matrix; } - EIGEN_DEVICE_FUNC const MatrixTypeNestedCleaned& nestedExpression() const { return m_matrix; } - EIGEN_DEVICE_FUNC MatrixTypeNestedCleaned& nestedExpression() { return m_matrix; } - - /** Efficient triangular matrix times vector/matrix product */ - template <typename OtherDerived> - EIGEN_DEVICE_FUNC const Product<SelfAdjointView, OtherDerived> operator*(const MatrixBase<OtherDerived>& rhs) const { - return Product<SelfAdjointView, OtherDerived>(*this, rhs.derived()); - } - - /** Efficient vector/matrix times triangular matrix product */ - template <typename OtherDerived> - friend EIGEN_DEVICE_FUNC const Product<OtherDerived, SelfAdjointView> operator*(const MatrixBase<OtherDerived>& lhs, - const SelfAdjointView& rhs) { - return Product<OtherDerived, SelfAdjointView>(lhs.derived(), rhs); - } + EIGEN_DEVICE_FUNC constexpr const MatrixTypeNestedCleaned& nestedExpression() const noexcept { return m_matrix; } + EIGEN_DEVICE_FUNC constexpr MatrixTypeNestedCleaned& nestedExpression() noexcept { return m_matrix; } EIGEN_DEVICE_FUNC const SelfAdjointView<const EIGEN_EXPR_BINARYOP_SCALAR_RETURN_TYPE(MatrixType, Scalar, product), UpLo> @@ -180,40 +190,6 @@ TriangularView<typename MatrixType::AdjointReturnType, TriMode> >(tmp2); } - typedef SelfAdjointView<const MatrixConjugateReturnType, UpLo> ConjugateReturnType; - /** \sa MatrixBase::conjugate() const */ - EIGEN_DEVICE_FUNC 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 std::conditional_t<Cond, ConjugateReturnType, ConstSelfAdjointView> conjugateIf() const { - typedef std::conditional_t<Cond, ConjugateReturnType, ConstSelfAdjointView> ReturnType; - return ReturnType(m_matrix.template conjugateIf<Cond>()); - } - - typedef SelfAdjointView<const typename MatrixType::AdjointReturnType, TransposeMode> AdjointReturnType; - /** \sa MatrixBase::adjoint() const */ - EIGEN_DEVICE_FUNC inline const AdjointReturnType adjoint() const { return AdjointReturnType(m_matrix.adjoint()); } - - typedef SelfAdjointView<typename MatrixType::TransposeReturnType, TransposeMode> TransposeReturnType; - /** \sa MatrixBase::transpose() */ - template <class Dummy = int> - EIGEN_DEVICE_FUNC inline TransposeReturnType transpose( - std::enable_if_t<Eigen::internal::is_lvalue<MatrixType>::value, Dummy*> = nullptr) { - typename MatrixType::TransposeReturnType tmp(m_matrix); - return TransposeReturnType(tmp); - } - - typedef SelfAdjointView<const typename MatrixType::ConstTransposeReturnType, TransposeMode> ConstTransposeReturnType; - /** \sa MatrixBase::transpose() const */ - EIGEN_DEVICE_FUNC inline const ConstTransposeReturnType transpose() const { - return ConstTransposeReturnType(m_matrix.transpose()); - } - /** \returns a const expression of the main diagonal of the matrix \c *this * * This method simply returns the diagonal of the nested expression, thus by-passing the SelfAdjointView decorator.
diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index 2b0f56b..fec3a2d 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h
@@ -21,8 +21,43 @@ template <int Side, typename TriangularType, typename Rhs> struct triangular_solve_retval; +template <typename MatrixType, unsigned int Mode, bool IsSelfAdjoint = (int(Mode) & int(SelfAdjoint)) != 0> +struct triangular_base_return_types { + typedef internal::remove_all_t<typename MatrixType::ConjugateReturnType> MatrixConjugateReturnType; + enum { + TransposeMode = (int(Mode) & int(Upper) ? Lower : 0) | (int(Mode) & int(Lower) ? Upper : 0) | + (int(Mode) & int(UnitDiag)) | (int(Mode) & int(ZeroDiag)) + }; + typedef TriangularView<std::add_const_t<MatrixType>, Mode> ConstView; + typedef TriangularView<const MatrixConjugateReturnType, Mode> ConjugateReturnType; + typedef TriangularView<const typename MatrixType::AdjointReturnType, TransposeMode> AdjointReturnType; + typedef TriangularView<typename MatrixType::TransposeReturnType, TransposeMode> TransposeReturnType; + typedef TriangularView<const typename MatrixType::ConstTransposeReturnType, TransposeMode> ConstTransposeReturnType; +}; + +template <typename MatrixType, unsigned int Mode> +struct triangular_base_return_types<MatrixType, Mode, true> { + typedef internal::remove_all_t<typename MatrixType::ConjugateReturnType> MatrixConjugateReturnType; + enum { + TriangularPart = int(Mode) & int(Upper | Lower), + TransposeMode = (int(Mode) & int(Upper) ? Lower : 0) | (int(Mode) & int(Lower) ? Upper : 0) + }; + typedef SelfAdjointView<std::add_const_t<MatrixType>, TriangularPart> ConstView; + typedef SelfAdjointView<const MatrixConjugateReturnType, TriangularPart> ConjugateReturnType; + typedef SelfAdjointView<const typename MatrixType::AdjointReturnType, TransposeMode> AdjointReturnType; + typedef SelfAdjointView<typename MatrixType::TransposeReturnType, TransposeMode> TransposeReturnType; + typedef SelfAdjointView<const typename MatrixType::ConstTransposeReturnType, TransposeMode> ConstTransposeReturnType; +}; + +template <unsigned int Mode, typename Expression> +EIGEN_DEVICE_FUNC inline typename triangular_base_return_types<Expression, Mode>::ConstView +make_triangular_base_cwise_view(const Expression& expression) { + typedef typename triangular_base_return_types<Expression, Mode>::ConstView ReturnType; + return ReturnType(expression); } +} // namespace internal + /** \class TriangularBase * \ingroup Core_Module * @@ -33,6 +68,8 @@ public: enum { Mode = internal::traits<Derived>::Mode, + TransposeMode = (int(Mode) & int(Upper) ? Lower : 0) | (int(Mode) & int(Lower) ? Upper : 0) | + (int(Mode) & int(UnitDiag)) | (int(Mode) & int(ZeroDiag)), RowsAtCompileTime = internal::traits<Derived>::RowsAtCompileTime, ColsAtCompileTime = internal::traits<Derived>::ColsAtCompileTime, MaxRowsAtCompileTime = internal::traits<Derived>::MaxRowsAtCompileTime, @@ -51,17 +88,25 @@ typedef typename internal::traits<Derived>::StorageKind StorageKind; typedef typename internal::traits<Derived>::StorageIndex StorageIndex; typedef typename internal::traits<Derived>::FullMatrixType DenseMatrixType; + typedef typename internal::traits<Derived>::ExpressionType ExpressionType; typedef DenseMatrixType DenseType; typedef Derived const& Nested; + typedef internal::triangular_base_return_types<ExpressionType, Mode> ReturnTypes; + typedef typename ReturnTypes::ConstView ConstView; + typedef typename ReturnTypes::ConjugateReturnType ConjugateReturnType; + typedef typename ReturnTypes::AdjointReturnType AdjointReturnType; + typedef typename ReturnTypes::TransposeReturnType TransposeReturnType; + typedef typename ReturnTypes::ConstTransposeReturnType ConstTransposeReturnType; EIGEN_DEVICE_FUNC inline TriangularBase() { eigen_assert(!((int(Mode) & int(UnitDiag)) && (int(Mode) & int(ZeroDiag)))); } - EIGEN_DEVICE_FUNC constexpr Index rows() const noexcept { return derived().rows(); } - EIGEN_DEVICE_FUNC constexpr Index cols() const noexcept { return derived().cols(); } - EIGEN_DEVICE_FUNC constexpr Index outerStride() const noexcept { return derived().outerStride(); } - EIGEN_DEVICE_FUNC constexpr Index innerStride() const noexcept { return derived().innerStride(); } + EIGEN_DEVICE_FUNC constexpr Index rows() const noexcept { return derived().nestedExpression().rows(); } + EIGEN_DEVICE_FUNC constexpr Index cols() const noexcept { return derived().nestedExpression().cols(); } + EIGEN_DEVICE_FUNC constexpr Index size() const noexcept { return rows() * cols(); } + EIGEN_DEVICE_FUNC constexpr Index outerStride() const noexcept { return derived().nestedExpression().outerStride(); } + EIGEN_DEVICE_FUNC constexpr Index innerStride() const noexcept { return derived().nestedExpression().innerStride(); } // dummy resize function EIGEN_DEVICE_FUNC void resize(Index rows, Index cols) { @@ -70,14 +115,48 @@ eigen_assert(rows == this->rows() && cols == this->cols()); } - EIGEN_DEVICE_FUNC inline Scalar coeff(Index row, Index col) const { return derived().coeff(row, col); } - EIGEN_DEVICE_FUNC inline Scalar& coeffRef(Index row, Index col) { return derived().coeffRef(row, col); } + /** \sa MatrixBase::fill() */ + EIGEN_DEVICE_FUNC void fill(const Scalar& value) { setConstant(value); } + + /** \sa MatrixBase::setConstant() */ + EIGEN_DEVICE_FUNC Derived& setConstant(const Scalar& value) { + EIGEN_STATIC_ASSERT_LVALUE(Derived); + return derived() = DenseMatrixType::Constant(rows(), cols(), value); + } + + /** \sa MatrixBase::setZero() */ + EIGEN_DEVICE_FUNC Derived& setZero() { return setConstant(Scalar(0)); } + + /** \sa MatrixBase::setOnes() */ + EIGEN_DEVICE_FUNC Derived& setOnes() { return setConstant(Scalar(1)); } + + /** \sa MatrixBase::setRandom() */ + EIGEN_DEVICE_FUNC Derived& setRandom() { + EIGEN_STATIC_ASSERT_LVALUE(Derived); + return derived() = DenseMatrixType::Random(rows(), cols()); + } + + /** \sa MatrixBase::setIdentity() */ + EIGEN_DEVICE_FUNC Derived& setIdentity() { + EIGEN_STATIC_ASSERT_LVALUE(Derived); + return derived() = DenseMatrixType::Identity(rows(), cols()); + } + + EIGEN_DEVICE_FUNC inline Scalar coeff(Index row, Index col) const { + check_coordinates_internal(row, col); + return derived().nestedExpression().coeff(row, col); + } + EIGEN_DEVICE_FUNC inline Scalar& coeffRef(Index row, Index col) { + EIGEN_STATIC_ASSERT_LVALUE(Derived); + check_coordinates_internal(row, col); + return derived().nestedExpression().coeffRef(row, col); + } /** \see MatrixBase::copyCoeff(row,col) */ template <typename Other> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void copyCoeff(Index row, Index col, Other& other) { - derived().coeffRef(row, col) = other.coeff(row, col); + coeffRef(row, col) = other.coeff(row, col); } EIGEN_DEVICE_FUNC inline Scalar operator()(Index row, Index col) const { @@ -95,8 +174,10 @@ #endif #ifndef EIGEN_PARSED_BY_DOXYGEN - EIGEN_DEVICE_FUNC inline const Derived& derived() const { return *static_cast<const Derived*>(this); } - EIGEN_DEVICE_FUNC inline Derived& derived() { return *static_cast<Derived*>(this); } + EIGEN_DEVICE_FUNC constexpr inline const Derived& derived() const noexcept { + return *static_cast<const Derived*>(this); + } + EIGEN_DEVICE_FUNC constexpr inline Derived& derived() noexcept { return *static_cast<Derived*>(this); } #endif // not EIGEN_PARSED_BY_DOXYGEN template <typename DenseDerived> @@ -104,6 +185,80 @@ template <typename DenseDerived> EIGEN_DEVICE_FUNC void evalToLazy(MatrixBase<DenseDerived>& other) const; + template <typename OtherDerived> + EIGEN_DEVICE_FUNC const Product<Derived, OtherDerived> operator*(const MatrixBase<OtherDerived>& rhs) const { + return Product<Derived, OtherDerived>(derived(), rhs.derived()); + } + + template <typename OtherDerived> + EIGEN_DEVICE_FUNC const Product<Derived, OtherDerived> operator*(const DiagonalBase<OtherDerived>& rhs) const { + return Product<Derived, OtherDerived>(derived(), rhs.derived()); + } + + template <typename OtherDerived, + std::enable_if_t<int(Mode) == int(OtherDerived::Mode) && (int(Mode) & int(UnitDiag)) == 0, int> = 0> + EIGEN_DEVICE_FUNC inline auto operator+(const TriangularBase<OtherDerived>& other) const + -> decltype(internal::make_triangular_base_cwise_view<Mode>( + std::declval<const ExpressionType&>() + + std::declval<const typename internal::traits<OtherDerived>::ExpressionType&>())) { + return internal::make_triangular_base_cwise_view<Mode>(derived().nestedExpression() + + other.derived().nestedExpression()); + } + + template <typename OtherDerived, + std::enable_if_t<int(Mode) == int(OtherDerived::Mode) && (int(Mode) & int(UnitDiag)) == 0, int> = 0> + EIGEN_DEVICE_FUNC inline auto operator-(const TriangularBase<OtherDerived>& other) const + -> decltype(internal::make_triangular_base_cwise_view<Mode>( + std::declval<const ExpressionType&>() - + std::declval<const typename internal::traits<OtherDerived>::ExpressionType&>())) { + return internal::make_triangular_base_cwise_view<Mode>(derived().nestedExpression() - + other.derived().nestedExpression()); + } + + template <typename OtherDerived> + friend EIGEN_DEVICE_FUNC const Product<OtherDerived, Derived> operator*(const MatrixBase<OtherDerived>& lhs, + const Derived& rhs) { + return Product<OtherDerived, Derived>(lhs.derived(), rhs); + } + + template <typename OtherDerived> + friend EIGEN_DEVICE_FUNC const Product<OtherDerived, Derived> operator*(const DiagonalBase<OtherDerived>& lhs, + const Derived& rhs) { + return Product<OtherDerived, Derived>(lhs.derived(), rhs); + } + + /** \sa MatrixBase::conjugate() const */ + EIGEN_DEVICE_FUNC inline const ConjugateReturnType conjugate() const { + return ConjugateReturnType(derived().nestedExpression().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 std::conditional_t<Cond, ConjugateReturnType, ConstView> conjugateIf() const { + typedef std::conditional_t<Cond, ConjugateReturnType, ConstView> ReturnType; + return ReturnType(derived().nestedExpression().template conjugateIf<Cond>()); + } + + /** \sa MatrixBase::adjoint() const */ + EIGEN_DEVICE_FUNC inline const AdjointReturnType adjoint() const { + return AdjointReturnType(derived().nestedExpression().adjoint()); + } + + /** \sa MatrixBase::transpose() */ + template <class Dummy = int> + EIGEN_DEVICE_FUNC inline TransposeReturnType transpose( + std::enable_if_t<Eigen::internal::is_lvalue<ExpressionType>::value, Dummy*> = nullptr) { + typename ExpressionType::TransposeReturnType tmp(derived().nestedExpression()); + return TransposeReturnType(tmp); + } + + /** \sa MatrixBase::transpose() const */ + EIGEN_DEVICE_FUNC inline const ConstTransposeReturnType transpose() const { + return ConstTransposeReturnType(derived().nestedExpression().transpose()); + } + EIGEN_DEVICE_FUNC DenseMatrixType toDenseMatrix() const { DenseMatrixType res(rows(), cols()); evalToLazy(res); @@ -179,9 +334,6 @@ typedef typename internal::traits<TriangularView>::MatrixTypeNested MatrixTypeNested; typedef typename internal::traits<TriangularView>::MatrixTypeNestedNonRef MatrixTypeNestedNonRef; - typedef internal::remove_all_t<typename MatrixType::ConjugateReturnType> MatrixConjugateReturnType; - typedef TriangularView<std::add_const_t<MatrixType>, Mode_> ConstTriangularView; - public: typedef typename internal::traits<TriangularView>::StorageKind StorageKind; typedef typename internal::traits<TriangularView>::MatrixTypeNestedCleaned NestedExpression; @@ -198,50 +350,11 @@ EIGEN_INHERIT_ASSIGNMENT_OPERATORS(TriangularView) - /** \copydoc EigenBase::rows() */ - EIGEN_DEVICE_FUNC constexpr Index rows() const noexcept { return m_matrix.rows(); } - /** \copydoc EigenBase::cols() */ - EIGEN_DEVICE_FUNC constexpr Index cols() const noexcept { return m_matrix.cols(); } - /** \returns a const reference to the nested expression */ - EIGEN_DEVICE_FUNC const NestedExpression& nestedExpression() const { return m_matrix; } + EIGEN_DEVICE_FUNC constexpr const NestedExpression& nestedExpression() const noexcept { return m_matrix; } /** \returns a reference to the nested expression */ - EIGEN_DEVICE_FUNC NestedExpression& nestedExpression() { return m_matrix; } - - typedef TriangularView<const MatrixConjugateReturnType, Mode> ConjugateReturnType; - /** \sa MatrixBase::conjugate() const */ - EIGEN_DEVICE_FUNC 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 std::conditional_t<Cond, ConjugateReturnType, ConstTriangularView> conjugateIf() const { - typedef std::conditional_t<Cond, ConjugateReturnType, ConstTriangularView> ReturnType; - return ReturnType(m_matrix.template conjugateIf<Cond>()); - } - - typedef TriangularView<const typename MatrixType::AdjointReturnType, TransposeMode> AdjointReturnType; - /** \sa MatrixBase::adjoint() const */ - EIGEN_DEVICE_FUNC inline const AdjointReturnType adjoint() const { return AdjointReturnType(m_matrix.adjoint()); } - - typedef TriangularView<typename MatrixType::TransposeReturnType, TransposeMode> TransposeReturnType; - /** \sa MatrixBase::transpose() */ - template <class Dummy = int> - EIGEN_DEVICE_FUNC inline TransposeReturnType transpose( - std::enable_if_t<Eigen::internal::is_lvalue<MatrixType>::value, Dummy*> = nullptr) { - typename MatrixType::TransposeReturnType tmp(m_matrix); - return TransposeReturnType(tmp); - } - - typedef TriangularView<const typename MatrixType::ConstTransposeReturnType, TransposeMode> ConstTransposeReturnType; - /** \sa MatrixBase::transpose() const */ - EIGEN_DEVICE_FUNC inline const ConstTransposeReturnType transpose() const { - return ConstTransposeReturnType(m_matrix.transpose()); - } + EIGEN_DEVICE_FUNC constexpr NestedExpression& nestedExpression() noexcept { return m_matrix; } template <typename Other> EIGEN_DEVICE_FUNC inline Solve<TriangularView, Other> solve(const MatrixBase<Other>& other) const { @@ -313,18 +426,12 @@ public: using Base::derived; using Base::evalToLazy; + using Base::operator*; typedef typename internal::traits<TriangularViewType>::StorageKind StorageKind; enum { Mode = Mode_, Flags = internal::traits<TriangularViewType>::Flags }; - /** \returns the outer-stride of the underlying dense matrix - * \sa DenseCoeffsBase::outerStride() */ - EIGEN_DEVICE_FUNC inline Index outerStride() const { return derived().nestedExpression().outerStride(); } - /** \returns the inner-stride of the underlying dense matrix - * \sa DenseCoeffsBase::innerStride() */ - EIGEN_DEVICE_FUNC inline Index innerStride() const { return derived().nestedExpression().innerStride(); } - /** \sa MatrixBase::operator+=() */ template <typename Other> EIGEN_DEVICE_FUNC TriangularViewType& operator+=(const DenseBase<Other>& other) { @@ -349,34 +456,6 @@ return *this = derived().nestedExpression() / other; } - /** \sa MatrixBase::fill() */ - EIGEN_DEVICE_FUNC void fill(const Scalar& value) { setConstant(value); } - /** \sa MatrixBase::setConstant() */ - EIGEN_DEVICE_FUNC TriangularViewType& setConstant(const Scalar& value) { - return *this = MatrixType::Constant(derived().rows(), derived().cols(), value); - } - /** \sa MatrixBase::setZero() */ - EIGEN_DEVICE_FUNC TriangularViewType& setZero() { return setConstant(Scalar(0)); } - /** \sa MatrixBase::setOnes() */ - EIGEN_DEVICE_FUNC TriangularViewType& setOnes() { return setConstant(Scalar(1)); } - - /** \sa MatrixBase::coeff() - * \warning the coordinates must fit into the referenced triangular part - */ - EIGEN_DEVICE_FUNC inline Scalar coeff(Index row, Index col) const { - Base::check_coordinates_internal(row, col); - return derived().nestedExpression().coeff(row, col); - } - - /** \sa MatrixBase::coeffRef() - * \warning the coordinates must fit into the referenced triangular part - */ - EIGEN_DEVICE_FUNC inline Scalar& coeffRef(Index row, Index col) { - EIGEN_STATIC_ASSERT_LVALUE(TriangularViewType); - Base::check_coordinates_internal(row, col); - return derived().nestedExpression().coeffRef(row, col); - } - /** Assigns a triangular matrix to a triangular part of a dense matrix */ template <typename OtherDerived> EIGEN_DEVICE_FUNC TriangularViewType& operator=(const TriangularBase<OtherDerived>& other); @@ -399,20 +478,6 @@ EIGEN_DEPRECATED EIGEN_DEVICE_FUNC void lazyAssign(const MatrixBase<OtherDerived>& other); #endif - /** Efficient triangular matrix times vector/matrix product */ - template <typename OtherDerived> - EIGEN_DEVICE_FUNC const Product<TriangularViewType, OtherDerived> operator*( - const MatrixBase<OtherDerived>& rhs) const { - return Product<TriangularViewType, OtherDerived>(derived(), rhs.derived()); - } - - /** Efficient vector/matrix times triangular matrix product */ - template <typename OtherDerived> - friend EIGEN_DEVICE_FUNC const Product<OtherDerived, TriangularViewType> operator*( - const MatrixBase<OtherDerived>& lhs, const TriangularViewImpl& rhs) { - return Product<OtherDerived, TriangularViewType>(lhs.derived(), rhs.derived()); - } - // Scaling a unit triangular view would break its implicit unit diagonal, so only non-unit modes participate. template <unsigned int M = Mode, std::enable_if_t<(M & UnitDiag) == 0, int> = 0> EIGEN_DEVICE_FUNC const @@ -772,6 +837,17 @@ typedef Dense2Triangular Kind; }; +template <typename Shape> +struct is_dense_structured_shape : std::integral_constant<bool, std::is_same<Shape, TriangularShape>::value || + std::is_same<Shape, SelfAdjointShape>::value> {}; + +template <typename Lhs, typename Rhs> +struct is_dense_structured_diagonal_product + : std::integral_constant<bool, (is_dense_structured_shape<typename evaluator_traits<Lhs>::Shape>::value && + std::is_same<typename evaluator_traits<Rhs>::Shape, DiagonalShape>::value) || + (std::is_same<typename evaluator_traits<Lhs>::Shape, DiagonalShape>::value && + is_dense_structured_shape<typename evaluator_traits<Rhs>::Shape>::value)> {}; + template <typename DstXprType, typename SrcXprType, typename Functor> struct Assignment<DstXprType, SrcXprType, Functor, Triangular2Triangular> { EIGEN_DEVICE_FUNC static void run(DstXprType& dst, const SrcXprType& src, const Functor& func) { @@ -897,18 +973,42 @@ namespace internal { +template <bool UseTriangularAssignmentLoop> +struct triangular_product_assignment_dispatcher { + template <typename DstXprType, typename SrcXprType, typename Functor, typename Scalar> + static void run(DstXprType& dst, const SrcXprType& src, const Functor&, const Scalar& alpha, bool beta) { + if (!beta) { + Index dstRows = src.rows(); + Index dstCols = src.cols(); + if ((dst.rows() != dstRows) || (dst.cols() != dstCols)) dst.resize(dstRows, dstCols); + } + + dst._assignProduct(src, alpha, beta); + } +}; + +template <> +struct triangular_product_assignment_dispatcher<true> { + template <typename DstXprType, typename SrcXprType, typename Functor, typename Scalar> + static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(DstXprType& dst, const SrcXprType& src, const Functor& func, + const Scalar& alpha, bool beta) { + EIGEN_UNUSED_VARIABLE(alpha); + EIGEN_UNUSED_VARIABLE(beta); + EIGEN_STATIC_ASSERT((int(DstXprType::Mode) & int(UnitDiag)) == 0, + WRITING_TO_TRIANGULAR_PART_WITH_UNIT_DIAGONAL_IS_NOT_SUPPORTED); + call_triangular_assignment_loop<DstXprType::Mode, false>(dst, src, func); + } +}; + // Triangular = Product template <typename DstXprType, typename Lhs, typename Rhs, typename Scalar> struct Assignment<DstXprType, Product<Lhs, Rhs, DefaultProduct>, internal::assign_op<Scalar, typename Product<Lhs, Rhs, DefaultProduct>::Scalar>, Dense2Triangular> { typedef Product<Lhs, Rhs, DefaultProduct> SrcXprType; static void run(DstXprType& dst, const SrcXprType& src, - const internal::assign_op<Scalar, typename SrcXprType::Scalar>&) { - Index dstRows = src.rows(); - Index dstCols = src.cols(); - if ((dst.rows() != dstRows) || (dst.cols() != dstCols)) dst.resize(dstRows, dstCols); - - dst._assignProduct(src, Scalar(1), false); + const internal::assign_op<Scalar, typename SrcXprType::Scalar>& func) { + enum { UseTriangularAssignmentLoop = is_dense_structured_diagonal_product<Lhs, Rhs>::value }; + triangular_product_assignment_dispatcher<UseTriangularAssignmentLoop>::run(dst, src, func, Scalar(1), false); } }; @@ -919,8 +1019,9 @@ Dense2Triangular> { typedef Product<Lhs, Rhs, DefaultProduct> SrcXprType; static void run(DstXprType& dst, const SrcXprType& src, - const internal::add_assign_op<Scalar, typename SrcXprType::Scalar>&) { - dst._assignProduct(src, Scalar(1), true); + const internal::add_assign_op<Scalar, typename SrcXprType::Scalar>& func) { + enum { UseTriangularAssignmentLoop = is_dense_structured_diagonal_product<Lhs, Rhs>::value }; + triangular_product_assignment_dispatcher<UseTriangularAssignmentLoop>::run(dst, src, func, Scalar(1), true); } }; @@ -931,8 +1032,9 @@ Dense2Triangular> { typedef Product<Lhs, Rhs, DefaultProduct> SrcXprType; static void run(DstXprType& dst, const SrcXprType& src, - const internal::sub_assign_op<Scalar, typename SrcXprType::Scalar>&) { - dst._assignProduct(src, Scalar(-1), true); + const internal::sub_assign_op<Scalar, typename SrcXprType::Scalar>& func) { + enum { UseTriangularAssignmentLoop = is_dense_structured_diagonal_product<Lhs, Rhs>::value }; + triangular_product_assignment_dispatcher<UseTriangularAssignmentLoop>::run(dst, src, func, Scalar(-1), true); } };
diff --git a/Eigen/src/SparseCore/SparseDiagonalProduct.h b/Eigen/src/SparseCore/SparseDiagonalProduct.h index e02fdaa..ae6f054 100644 --- a/Eigen/src/SparseCore/SparseDiagonalProduct.h +++ b/Eigen/src/SparseCore/SparseDiagonalProduct.h
@@ -48,7 +48,7 @@ typedef sparse_diagonal_product_evaluator<Rhs, typename Lhs::DiagonalVectorType, Rhs::Flags & RowMajorBit ? SDP_AsScalarProduct : SDP_AsCwiseProduct> Base; - explicit product_evaluator(const XprType &xpr) : Base(xpr.rhs(), xpr.lhs().diagonal()) {} + explicit product_evaluator(const XprType& xpr) : Base(xpr.rhs(), xpr.lhs().diagonal()) {} }; template <typename Lhs, typename Rhs, int ProductTag> @@ -65,7 +65,149 @@ typedef sparse_diagonal_product_evaluator<Lhs, Transpose<const typename Rhs::DiagonalVectorType>, Lhs::Flags & RowMajorBit ? SDP_AsCwiseProduct : SDP_AsScalarProduct> Base; - explicit product_evaluator(const XprType &xpr) : Base(xpr.lhs(), xpr.rhs().diagonal().transpose()) {} + explicit product_evaluator(const XprType& xpr) : Base(xpr.lhs(), xpr.rhs().diagonal().transpose()) {} +}; + +// SparseSelfAdjointView synthesizes mirrored entries; build the diagonal-scaled +// full result directly instead of first materializing an unscaled PlainObject. +template <int Mode, int ProductOrder, typename SelfAdjointViewType, typename DiagonalType, typename Dest> +struct sparse_selfadjoint_diagonal_product_impl { + typedef typename SelfAdjointViewType::MatrixTypeNested_ MatrixType; + typedef evaluator<MatrixType> MatrixEvaluator; + typedef typename MatrixEvaluator::InnerIterator MatrixIterator; + typedef typename Dest::StorageIndex StorageIndex; + typedef Matrix<StorageIndex, Dynamic, 1> VectorI; + enum { IsFullMode = Mode == int(Upper | Lower), IsLowerMode = (Mode & int(Lower)) == int(Lower) }; + + static void run(Dest& dest, const SelfAdjointViewType& selfadjoint, const DiagonalType& diagonal) { + MatrixEvaluator matrixEval(selfadjoint.matrix()); + const Index size = selfadjoint.rows(); + VectorI count(size); + count.setZero(); + dest.resize(size, size); + + for (Index outer = 0; outer < selfadjoint.matrix().outerSize(); ++outer) { + for (MatrixIterator it(matrixEval, outer); it; ++it) countEntry(count, it.row(), it.col()); + } + + Index nnz = count.sum(); + dest.resizeNonZeros(nnz); + dest.outerIndexPtr()[0] = 0; + for (Index outer = 0; outer < size; ++outer) + dest.outerIndexPtr()[outer + 1] = dest.outerIndexPtr()[outer] + count[outer]; + for (Index outer = 0; outer < size; ++outer) count[outer] = dest.outerIndexPtr()[outer]; + + for (Index outer = 0; outer < selfadjoint.matrix().outerSize(); ++outer) { + for (MatrixIterator it(matrixEval, outer); it; ++it) { + const Index row = it.row(); + const Index col = it.col(); + if (isStored(row, col)) { + insertEntry(dest, count, diagonal, row, col, it.value()); + } + if (mirrorsStoredEntry(row, col)) { + insertEntry(dest, count, diagonal, col, row, numext::conj(it.value())); + } + } + } + } + + private: + static EIGEN_STRONG_INLINE bool isStored(Index row, Index col) { + return IsFullMode || row == col || (IsLowerMode ? row > col : row < col); + } + + static EIGEN_STRONG_INLINE bool mirrorsStoredEntry(Index row, Index col) { + return !IsFullMode && row != col && (IsLowerMode ? row > col : row < col); + } + + static void countEntry(VectorI& count, Index row, Index col) { + if (isStored(row, col)) { + ++count[outerIndex(row, col)]; + } + if (mirrorsStoredEntry(row, col)) { + ++count[outerIndex(col, row)]; + } + } + + static EIGEN_STRONG_INLINE StorageIndex outerIndex(Index row, Index col) { + return internal::convert_index<StorageIndex>(Dest::IsRowMajor ? row : col); + } + + static EIGEN_STRONG_INLINE StorageIndex innerIndex(Index row, Index col) { + return internal::convert_index<StorageIndex>(Dest::IsRowMajor ? col : row); + } + + template <typename Coeff> + static void insertEntry(Dest& dest, VectorI& count, const DiagonalType& diagonal, Index row, Index col, + const Coeff& coeff) { + const StorageIndex outer = outerIndex(row, col); + const Index k = count[outer]++; + dest.innerIndexPtr()[k] = innerIndex(row, col); + if (ProductOrder == OnTheLeft) + dest.valuePtr()[k] = diagonal.coeff(row) * coeff; + else + dest.valuePtr()[k] = coeff * diagonal.coeff(col); + } +}; + +template <typename Lhs, typename Rhs> +struct materialized_left_sparse_product_evaluator_base + : public evaluator<typename Product<Lhs, typename Rhs::PlainObject, DefaultProduct>::PlainObject> { + typedef Product<Lhs, Rhs, DefaultProduct> XprType; + typedef typename XprType::PlainObject PlainObject; + typedef evaluator<PlainObject> Base; + + explicit materialized_left_sparse_product_evaluator_base(const XprType& xpr) : m_result(xpr.rows(), xpr.cols()) { + internal::construct_at<Base>(this, m_result); + sparse_selfadjoint_diagonal_product_impl<Rhs::Mode, OnTheLeft, Rhs, typename Lhs::DiagonalVectorType, + PlainObject>::run(m_result, xpr.rhs(), xpr.lhs().diagonal()); + } + + protected: + PlainObject m_result; +}; + +template <typename Lhs, typename Rhs> +struct materialized_right_sparse_product_evaluator_base + : public evaluator<typename Product<typename Lhs::PlainObject, Rhs, DefaultProduct>::PlainObject> { + typedef Product<Lhs, Rhs, DefaultProduct> XprType; + typedef typename XprType::PlainObject PlainObject; + typedef evaluator<PlainObject> Base; + + explicit materialized_right_sparse_product_evaluator_base(const XprType& xpr) : m_result(xpr.rows(), xpr.cols()) { + internal::construct_at<Base>(this, m_result); + sparse_selfadjoint_diagonal_product_impl<Lhs::Mode, OnTheRight, Lhs, typename Rhs::DiagonalVectorType, + PlainObject>::run(m_result, xpr.lhs(), xpr.rhs().diagonal()); + } + + protected: + PlainObject m_result; +}; + +template <typename Lhs, typename Rhs, int ProductTag> +struct product_evaluator<Product<Lhs, Rhs, DefaultProduct>, ProductTag, DiagonalShape, SparseTriangularShape> + : product_evaluator<Product<Lhs, Rhs, DefaultProduct>, ProductTag, DiagonalShape, SparseShape> { + typedef product_evaluator<Product<Lhs, Rhs, DefaultProduct>, ProductTag, DiagonalShape, SparseShape> Base; + using Base::Base; +}; + +template <typename Lhs, typename Rhs, int ProductTag> +struct product_evaluator<Product<Lhs, Rhs, DefaultProduct>, ProductTag, DiagonalShape, SparseSelfAdjointShape> + : materialized_left_sparse_product_evaluator_base<Lhs, Rhs> { + using materialized_left_sparse_product_evaluator_base<Lhs, Rhs>::materialized_left_sparse_product_evaluator_base; +}; + +template <typename Lhs, typename Rhs, int ProductTag> +struct product_evaluator<Product<Lhs, Rhs, DefaultProduct>, ProductTag, SparseTriangularShape, DiagonalShape> + : product_evaluator<Product<Lhs, Rhs, DefaultProduct>, ProductTag, SparseShape, DiagonalShape> { + typedef product_evaluator<Product<Lhs, Rhs, DefaultProduct>, ProductTag, SparseShape, DiagonalShape> Base; + using Base::Base; +}; + +template <typename Lhs, typename Rhs, int ProductTag> +struct product_evaluator<Product<Lhs, Rhs, DefaultProduct>, ProductTag, SparseSelfAdjointShape, DiagonalShape> + : materialized_right_sparse_product_evaluator_base<Lhs, Rhs> { + using materialized_right_sparse_product_evaluator_base<Lhs, Rhs>::materialized_right_sparse_product_evaluator_base; }; template <typename SparseXprType, typename DiagonalCoeffType> @@ -77,7 +219,7 @@ public: class InnerIterator : public SparseXprInnerIterator { public: - InnerIterator(const sparse_diagonal_product_evaluator &xprEval, Index outer) + InnerIterator(const sparse_diagonal_product_evaluator& xprEval, Index outer) : SparseXprInnerIterator(xprEval.m_sparseXprImpl, outer), m_coeff(xprEval.m_diagCoeffImpl.coeff(outer)) {} EIGEN_STRONG_INLINE Scalar value() const { return m_coeff * SparseXprInnerIterator::value(); } @@ -86,7 +228,7 @@ typename DiagonalCoeffType::Scalar m_coeff; }; - sparse_diagonal_product_evaluator(const SparseXprType &sparseXpr, const DiagonalCoeffType &diagCoeff) + sparse_diagonal_product_evaluator(const SparseXprType& sparseXpr, const DiagonalCoeffType& diagCoeff) : m_sparseXprImpl(sparseXpr), m_diagCoeffImpl(diagCoeff) {} Index nonZerosEstimate() const { return m_sparseXprImpl.nonZerosEstimate(); } @@ -109,7 +251,7 @@ typedef typename evaluator<SparseXprType>::InnerIterator SparseXprIter; public: - InnerIterator(const sparse_diagonal_product_evaluator &xprEval, Index outer) + InnerIterator(const sparse_diagonal_product_evaluator& xprEval, Index outer) : m_sparseIter(xprEval.m_sparseXprEval, outer), m_diagCoeffNested(xprEval.m_diagCoeffNested) {} inline Scalar value() const { return m_sparseIter.value() * m_diagCoeffNested.coeff(index()); } @@ -118,7 +260,7 @@ inline Index col() const { return SparseXprType::IsRowMajor ? m_sparseIter.index() : m_sparseIter.outer(); } inline Index row() const { return SparseXprType::IsRowMajor ? m_sparseIter.outer() : m_sparseIter.index(); } - EIGEN_STRONG_INLINE InnerIterator &operator++() { + EIGEN_STRONG_INLINE InnerIterator& operator++() { ++m_sparseIter; return *this; } @@ -129,7 +271,7 @@ DiagCoeffNested m_diagCoeffNested; }; - sparse_diagonal_product_evaluator(const SparseXprType &sparseXpr, const DiagCoeffType &diagCoeff) + sparse_diagonal_product_evaluator(const SparseXprType& sparseXpr, const DiagCoeffType& diagCoeff) : m_sparseXprEval(sparseXpr), m_diagCoeffNested(diagCoeff) {} Index nonZerosEstimate() const { return m_sparseXprEval.nonZerosEstimate(); }
diff --git a/Eigen/src/SparseCore/SparseSelfAdjointView.h b/Eigen/src/SparseCore/SparseSelfAdjointView.h index 32fbcd4..55326e9 100644 --- a/Eigen/src/SparseCore/SparseSelfAdjointView.h +++ b/Eigen/src/SparseCore/SparseSelfAdjointView.h
@@ -109,6 +109,11 @@ return Product<SparseSelfAdjointView, OtherDerived>(*this, rhs.derived()); } + template <typename OtherDerived> + Product<SparseSelfAdjointView, OtherDerived> operator*(const DiagonalBase<OtherDerived>& rhs) const { + return Product<SparseSelfAdjointView, OtherDerived>(*this, rhs.derived()); + } + /** Efficient dense vector/matrix times sparse self-adjoint matrix product */ template <typename OtherDerived> friend Product<OtherDerived, SparseSelfAdjointView> operator*(const MatrixBase<OtherDerived>& lhs, @@ -116,6 +121,12 @@ return Product<OtherDerived, SparseSelfAdjointView>(lhs.derived(), rhs); } + template <typename OtherDerived> + friend Product<OtherDerived, SparseSelfAdjointView> operator*(const DiagonalBase<OtherDerived>& lhs, + const SparseSelfAdjointView& rhs) { + return Product<OtherDerived, SparseSelfAdjointView>(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; }
diff --git a/Eigen/src/SparseCore/SparseTriangularView.h b/Eigen/src/SparseCore/SparseTriangularView.h index a6c3eaa..c98ef0f 100644 --- a/Eigen/src/SparseCore/SparseTriangularView.h +++ b/Eigen/src/SparseCore/SparseTriangularView.h
@@ -58,6 +58,9 @@ this->solveInPlace(dst); } + EIGEN_DEVICE_FUNC constexpr Index rows() const noexcept { return derived().nestedExpression().rows(); } + EIGEN_DEVICE_FUNC constexpr Index cols() const noexcept { return derived().nestedExpression().cols(); } + /** Applies the inverse of \c *this to the dense vector or matrix \a other, "in-place" */ template <typename OtherDerived> void solveInPlace(MatrixBase<OtherDerived>& other) const;
diff --git a/test/diagonalmatrices.cpp b/test/diagonalmatrices.cpp index 62450ab..e61b4b6 100644 --- a/test/diagonalmatrices.cpp +++ b/test/diagonalmatrices.cpp
@@ -172,6 +172,126 @@ VERIFY_IS_APPROX((res1 = points.topLeftCorner<2, 2>() * diag.asDiagonal()), res2 = tmp2 * diag.asDiagonal()); } +template <int> +void bug2013() { + Matrix3d m = Matrix3d::Random(); + Vector3d d = Vector3d::Random(); + + Matrix3d ref_unit_lower = m.template triangularView<UnitLower>(); + Matrix3d ref_unit_upper = m.template triangularView<UnitUpper>(); + Matrix3d ref_lower = m.template triangularView<Lower>(); + Matrix3d ref_upper = m.template triangularView<Upper>(); + + VERIFY_IS_APPROX((m.template triangularView<UnitLower>() * d.asDiagonal()).eval(), ref_unit_lower * d.asDiagonal()); + VERIFY_IS_APPROX((d.asDiagonal() * m.template triangularView<UnitLower>()).eval(), d.asDiagonal() * ref_unit_lower); + + VERIFY_IS_APPROX((m.template triangularView<UnitUpper>() * d.asDiagonal()).eval(), ref_unit_upper * d.asDiagonal()); + VERIFY_IS_APPROX((d.asDiagonal() * m.template triangularView<UnitUpper>()).eval(), d.asDiagonal() * ref_unit_upper); + + Matrix3d actual = Matrix3d::Random(); + Matrix3d expected = actual; + actual = m; + expected = m; + actual.template triangularView<Upper>() = actual.template triangularView<Upper>() * d.asDiagonal(); + expected.template triangularView<Upper>() = Matrix3d(ref_upper * d.asDiagonal()); + VERIFY_IS_APPROX(actual, expected); + + actual = m; + expected = m; + actual.template triangularView<Lower>() = d.asDiagonal() * actual.template triangularView<Lower>(); + expected.template triangularView<Lower>() = Matrix3d(d.asDiagonal() * ref_lower); + VERIFY_IS_APPROX(actual, expected); + + actual.setRandom(); + expected = actual; + actual.noalias() += m.template triangularView<UnitLower>() * d.asDiagonal(); + expected.noalias() += ref_unit_lower * d.asDiagonal(); + VERIFY_IS_APPROX(actual, expected); + + actual.setRandom(); + expected = actual; + actual.noalias() -= d.asDiagonal() * m.template triangularView<UnitUpper>(); + expected.noalias() -= d.asDiagonal() * ref_unit_upper; + VERIFY_IS_APPROX(actual, expected); + + MatrixXd dynamic_m = MatrixXd::Random(4, 4); + VectorXd dynamic_d = VectorXd::Random(4); + MatrixXd dynamic_expected(4, 4); + MatrixXd no_malloc_result(4, 4); + + dynamic_expected = MatrixXd(dynamic_m.template triangularView<UnitLower>()) * dynamic_d.asDiagonal(); + internal::set_is_malloc_allowed(false); + no_malloc_result.noalias() = dynamic_m.template triangularView<UnitLower>() * dynamic_d.asDiagonal(); + internal::set_is_malloc_allowed(true); + VERIFY_IS_APPROX(no_malloc_result, dynamic_expected); + + dynamic_expected = dynamic_d.asDiagonal() * MatrixXd(dynamic_m.template triangularView<UnitUpper>()); + internal::set_is_malloc_allowed(false); + 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); +} + +template <int> +void selfadjoint_diagonal_products() { + Matrix3cd m = Matrix3cd::Random(); + m.diagonal() = m.diagonal().real(); + Vector3cd d = Vector3cd::Random(); + + Matrix3cd ref_lower = m.template selfadjointView<Lower>(); + Matrix3cd ref_upper = m.template selfadjointView<Upper>(); + + VERIFY_IS_APPROX((m.template selfadjointView<Lower>() * d.asDiagonal()).eval(), ref_lower * d.asDiagonal()); + VERIFY_IS_APPROX((d.asDiagonal() * m.template selfadjointView<Lower>()).eval(), d.asDiagonal() * ref_lower); + + VERIFY_IS_APPROX((m.template selfadjointView<Upper>() * d.asDiagonal()).eval(), ref_upper * d.asDiagonal()); + VERIFY_IS_APPROX((d.asDiagonal() * m.template selfadjointView<Upper>()).eval(), d.asDiagonal() * ref_upper); + + Matrix3cd actual = Matrix3cd::Random(); + Matrix3cd expected = actual; + actual = m; + expected = m; + actual.template selfadjointView<Upper>() = actual.template selfadjointView<Upper>() * d.asDiagonal(); + expected.template triangularView<Upper>() = Matrix3cd(ref_upper * d.asDiagonal()); + VERIFY_IS_APPROX(actual, expected); + + actual = m; + expected = m; + actual.template selfadjointView<Lower>() = d.asDiagonal() * actual.template selfadjointView<Lower>(); + expected.template triangularView<Lower>() = Matrix3cd(d.asDiagonal() * ref_lower); + VERIFY_IS_APPROX(actual, expected); + + actual.setRandom(); + expected = actual; + actual.noalias() += m.template selfadjointView<Lower>() * d.asDiagonal(); + expected.noalias() += ref_lower * d.asDiagonal(); + VERIFY_IS_APPROX(actual, expected); + + actual.setRandom(); + expected = actual; + actual.noalias() -= d.asDiagonal() * m.template selfadjointView<Upper>(); + expected.noalias() -= d.asDiagonal() * ref_upper; + VERIFY_IS_APPROX(actual, expected); + + MatrixXcd dynamic_m = MatrixXcd::Random(4, 4); + dynamic_m.diagonal() = dynamic_m.diagonal().real(); + VectorXcd dynamic_d = VectorXcd::Random(4); + MatrixXcd dynamic_expected(4, 4); + MatrixXcd no_malloc_result(4, 4); + + dynamic_expected = MatrixXcd(dynamic_m.template selfadjointView<Lower>()) * dynamic_d.asDiagonal(); + internal::set_is_malloc_allowed(false); + no_malloc_result.noalias() = dynamic_m.template selfadjointView<Lower>() * dynamic_d.asDiagonal(); + internal::set_is_malloc_allowed(true); + VERIFY_IS_APPROX(no_malloc_result, dynamic_expected); + + dynamic_expected = dynamic_d.asDiagonal() * MatrixXcd(dynamic_m.template selfadjointView<Upper>()); + internal::set_is_malloc_allowed(false); + 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); +} + EIGEN_DECLARE_TEST(diagonalmatrices) { for (int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1(diagonalmatrices(Matrix<float, 1, 1>())); @@ -194,4 +314,6 @@ CALL_SUBTEST_9(as_scalar_product(MatrixXf(1, 1))); } CALL_SUBTEST_10(bug987<0>()); + CALL_SUBTEST_10(bug2013<0>()); + CALL_SUBTEST_10(selfadjoint_diagonal_products<0>()); }
diff --git a/test/product_extra.cpp b/test/product_extra.cpp index 4751b0d..9b6a330 100644 --- a/test/product_extra.cpp +++ b/test/product_extra.cpp
@@ -217,6 +217,13 @@ } template <int> +void triangular_product_assignment_size_mismatch() { + MatrixXd m1 = MatrixXd::Random(5, 5); + MatrixXd m2(6, 6); + VERIFY_RAISES_ASSERT(m2.triangularView<Lower>() = m1 * m1); +} + +template <int> void unaligned_objects() { // Regression test for the bug reported here: // http://forum.kde.org/viewtopic.php?f=74&t=107541 @@ -734,6 +741,7 @@ MatrixXf(internal::random<int>(1, EIGEN_TEST_MAX_SIZE), internal::random<int>(1, EIGEN_TEST_MAX_SIZE)))); } CALL_SUBTEST_5(bug_127<0>()); + CALL_SUBTEST_5(triangular_product_assignment_size_mismatch<0>()); CALL_SUBTEST_5(bug_817<0>()); CALL_SUBTEST_5(bug_1308<0>()); CALL_SUBTEST_6(unaligned_objects<0>());
diff --git a/test/product_notemporary.cpp b/test/product_notemporary.cpp index 27136cb..961d87a 100644 --- a/test/product_notemporary.cpp +++ b/test/product_notemporary.cpp
@@ -108,6 +108,11 @@ VERIFY_EVALUATION_COUNT(m3.template triangularView<Upper>() = (m1 * m2.adjoint()), 0); VERIFY_EVALUATION_COUNT(m3.template triangularView<Upper>() -= (m1 * m2.adjoint()), 0); + VERIFY_EVALUATION_COUNT(m3.template selfadjointView<Lower>() = m1.template triangularView<Upper>(), 0); + VERIFY_EVALUATION_COUNT(m3.template selfadjointView<Upper>() = m1.template selfadjointView<Lower>(), 0); + VERIFY_EVALUATION_COUNT(m3.template selfadjointView<Lower>() = + s1 * m1.template selfadjointView<Lower>() + m2.template selfadjointView<Lower>(), + 0); // NOTE this is because the blas_traits require innerstride==1 to avoid a temporary, but that doesn't seem to be // actually needed for the triangular products
diff --git a/test/product_selfadjoint.cpp b/test/product_selfadjoint.cpp index 4e45fe8..5eda21d 100644 --- a/test/product_selfadjoint.cpp +++ b/test/product_selfadjoint.cpp
@@ -29,6 +29,99 @@ m1 = (m1.adjoint() + m1).eval(); + // Dense selfadjoint assignment is documented as writing only the referenced triangle. + m2.setRandom(); + m3 = m2; + m2.template selfadjointView<Upper>() = m1 + m1.adjoint(); + m3.template triangularView<Upper>() = m1 + m1.adjoint(); + VERIFY_IS_APPROX(m2, m3); + + m2.setRandom(); + m3 = m2; + m2.template selfadjointView<Upper>().setZero(); + m3.template triangularView<Upper>().setZero(); + VERIFY_IS_APPROX(m2, m3); + + m2.setRandom(); + m3 = m2; + m2.template selfadjointView<Lower>().setOnes(); + m3.template triangularView<Lower>().setOnes(); + VERIFY_IS_APPROX(m2, m3); + + m2.setRandom(); + m3 = m2; + m2.template selfadjointView<Upper>().setConstant(s1); + m3.template triangularView<Upper>().setConstant(s1); + VERIFY_IS_APPROX(m2, m3); + + m2.setRandom(); + m3 = m2; + m2.template selfadjointView<Lower>().fill(s2); + m3.template triangularView<Lower>().fill(s2); + VERIFY_IS_APPROX(m2, m3); + + m2.setRandom(); + m3 = m2; + m2.template selfadjointView<Upper>().setIdentity(); + m3.template triangularView<Upper>().setIdentity(); + VERIFY_IS_APPROX(m2, m3); + + m2.setRandom(); + m3 = m2; + m2.template selfadjointView<Lower>().setRandom(); + VERIFY_IS_APPROX(m2.template triangularView<StrictlyUpper>().toDenseMatrix(), + m3.template triangularView<StrictlyUpper>().toDenseMatrix()); + + m2.setRandom(); + m3 = m2; + m2.template selfadjointView<Upper>() += m1; + m3.template triangularView<Upper>() += m1; + VERIFY_IS_APPROX(m2, m3); + + m2.setRandom(); + m3 = m2; + m2.template selfadjointView<Lower>() -= m1; + m3.template triangularView<Lower>() -= m1; + VERIFY_IS_APPROX(m2, m3); + + m2.setRandom(); + m3 = m2; + m2.template selfadjointView<Upper>() *= s1; + m3.template triangularView<Upper>() *= 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; + VERIFY_IS_APPROX(m2, m3); + + m2.setRandom(); + m2.template selfadjointView<Lower>() = m1.template selfadjointView<Lower>(); + m3 = m1.template selfadjointView<Lower>().toDenseMatrix(); + VERIFY_IS_APPROX(m2, m3); + + m2.setRandom(); + m2.template selfadjointView<Lower>() = m1.template selfadjointView<Upper>(); + m3 = m1.template selfadjointView<Upper>().toDenseMatrix(); + VERIFY_IS_APPROX(m2, m3); + + m2.setRandom(); + m2.template selfadjointView<Lower>() = m1.template triangularView<Upper>(); + m3 = m1.template triangularView<Upper>().toDenseMatrix(); + VERIFY_IS_APPROX(m2, m3); + + m2.setRandom(); + m2.template selfadjointView<Lower>() = m1.template triangularView<UnitUpper>(); + m3 = m1.template triangularView<UnitUpper>().toDenseMatrix(); + VERIFY_IS_APPROX(m2, m3); + + m2.setRandom(); + m2.template selfadjointView<Upper>() = m1.template triangularView<StrictlyLower>(); + m3 = m1.template triangularView<StrictlyLower>().toDenseMatrix(); + VERIFY_IS_APPROX(m2, m3); + // rank2 update m2 = m1.template triangularView<Lower>(); m2.template selfadjointView<Lower>().rankUpdate(v1, v2);
diff --git a/test/sparse_product.cpp b/test/sparse_product.cpp index 426a375..f001a90 100644 --- a/test/sparse_product.cpp +++ b/test/sparse_product.cpp
@@ -33,6 +33,37 @@ VERIFY((#XPR) && nb_temporaries == N); \ } +template <typename Lhs, typename Rhs, typename = void> +struct has_product : std::false_type {}; + +template <typename Lhs, typename Rhs> +struct has_product<Lhs, Rhs, internal::void_t<decltype(std::declval<const Lhs&>() * std::declval<const Rhs&>())>> + : std::true_type {}; + +template <typename SparseMatrixType> +void sparse_structured_view_product_sfinae() { + typedef typename SparseMatrixType::Scalar Scalar; + typedef Matrix<Scalar, Dynamic, Dynamic> DenseMatrixType; + typedef Matrix<Scalar, Dynamic, 1> DenseVectorType; + typedef decltype(std::declval<SparseMatrixType&>().template triangularView<Lower>()) TriangularViewType; + typedef decltype(std::declval<SparseMatrixType&>().template selfadjointView<Lower>()) SelfAdjointViewType; + typedef decltype(std::declval<const DenseVectorType&>().asDiagonal()) DiagonalType; + + STATIC_CHECK((has_product<TriangularViewType, SparseMatrixType>::value)); + STATIC_CHECK((has_product<SparseMatrixType, TriangularViewType>::value)); + STATIC_CHECK((has_product<TriangularViewType, DenseMatrixType>::value)); + STATIC_CHECK((has_product<DenseMatrixType, TriangularViewType>::value)); + STATIC_CHECK((has_product<TriangularViewType, DiagonalType>::value)); + STATIC_CHECK((has_product<DiagonalType, TriangularViewType>::value)); + + STATIC_CHECK((has_product<SelfAdjointViewType, SparseMatrixType>::value)); + STATIC_CHECK((has_product<SparseMatrixType, SelfAdjointViewType>::value)); + STATIC_CHECK((has_product<SelfAdjointViewType, DenseMatrixType>::value)); + STATIC_CHECK((has_product<DenseMatrixType, SelfAdjointViewType>::value)); + STATIC_CHECK((has_product<SelfAdjointViewType, DiagonalType>::value)); + STATIC_CHECK((has_product<DiagonalType, SelfAdjointViewType>::value)); +} + template <typename SparseMatrixType> void sparse_product() { typedef typename SparseMatrixType::StorageIndex StorageIndex; @@ -334,12 +365,26 @@ VERIFY_IS_APPROX(x.noalias() -= mLo.template selfadjointView<Lower>() * b, refX -= refS * b); VERIFY_IS_APPROX(x.noalias() += mS.template selfadjointView<Upper | Lower>() * b, refX += refS * b); + DenseVector scale = DenseVector::Random(rows); + VERIFY_IS_APPROX(x = mLo.template selfadjointView<Lower>() * scale.asDiagonal(), refX = refS * scale.asDiagonal()); + 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); + // sparse selfadjointView with sparse matrices SparseMatrixType mSres(rows, rows); VERIFY_IS_APPROX(mSres = mLo.template selfadjointView<Lower>() * mS, refX = refLo.template selfadjointView<Lower>() * refS); VERIFY_IS_APPROX(mSres = mS * mLo.template selfadjointView<Lower>(), refX = refS * refLo.template selfadjointView<Lower>()); + VERIFY_IS_APPROX(mSres = mLo.template selfadjointView<Lower>() * scale.asDiagonal(), + refX = refS * scale.asDiagonal()); + VERIFY_IS_APPROX(mSres = scale.asDiagonal() * mLo.template selfadjointView<Lower>(), + refX = scale.asDiagonal() * refS); + VERIFY_IS_APPROX(mSres = mUp.template selfadjointView<Upper>() * scale.asDiagonal(), + refX = refS * scale.asDiagonal()); + VERIFY_IS_APPROX(mSres = scale.asDiagonal() * mUp.template selfadjointView<Upper>(), + refX = scale.asDiagonal() * refS); // sparse triangularView with dense matrices VERIFY_IS_APPROX(x = mA.template triangularView<Upper>() * b, refX = refA.template triangularView<Upper>() * b); @@ -356,6 +401,15 @@ refX = refA.template triangularView<Upper>() * refS); VERIFY_IS_APPROX(mSres = mS * mA.template triangularView<Upper>(), refX = refS * refA.template triangularView<Upper>()); + + VERIFY_IS_APPROX(mSres = mA.template triangularView<UnitLower>() * scale.asDiagonal(), + refX = DenseMatrix(refA.template triangularView<UnitLower>()) * scale.asDiagonal()); + VERIFY_IS_APPROX(mSres = scale.asDiagonal() * mA.template triangularView<UnitLower>(), + refX = scale.asDiagonal() * DenseMatrix(refA.template triangularView<UnitLower>())); + VERIFY_IS_APPROX(mSres = mA.template triangularView<UnitUpper>() * scale.asDiagonal(), + refX = DenseMatrix(refA.template triangularView<UnitUpper>()) * scale.asDiagonal()); + VERIFY_IS_APPROX(mSres = scale.asDiagonal() * mA.template triangularView<UnitUpper>(), + refX = scale.asDiagonal() * DenseMatrix(refA.template triangularView<UnitUpper>())); } } @@ -554,16 +608,19 @@ } EIGEN_DECLARE_TEST(sparse_product) { + sparse_structured_view_product_sfinae<SparseMatrix<double, ColMajor>>(); + sparse_structured_view_product_sfinae<SparseMatrix<double, RowMajor>>(); + for (int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1((test_sparse_vector_dense_product())); - CALL_SUBTEST_1((sparse_product<SparseMatrix<double, ColMajor> >())); - CALL_SUBTEST_1((sparse_product<SparseMatrix<double, RowMajor> >())); + CALL_SUBTEST_1((sparse_product<SparseMatrix<double, ColMajor>>())); + CALL_SUBTEST_1((sparse_product<SparseMatrix<double, RowMajor>>())); CALL_SUBTEST_1((bug_942<double>())); - CALL_SUBTEST_2((sparse_product<SparseMatrix<std::complex<double>, ColMajor> >())); - CALL_SUBTEST_2((sparse_product<SparseMatrix<std::complex<double>, RowMajor> >())); - CALL_SUBTEST_3((sparse_product<SparseMatrix<float, ColMajor, long int> >())); - CALL_SUBTEST_4(( - sparse_product_regression_test<SparseMatrix<double, RowMajor>, Matrix<double, Dynamic, Dynamic, RowMajor> >())); + CALL_SUBTEST_2((sparse_product<SparseMatrix<std::complex<double>, ColMajor>>())); + CALL_SUBTEST_2((sparse_product<SparseMatrix<std::complex<double>, RowMajor>>())); + CALL_SUBTEST_3((sparse_product<SparseMatrix<float, ColMajor, long int>>())); + CALL_SUBTEST_4( + (sparse_product_regression_test<SparseMatrix<double, RowMajor>, Matrix<double, Dynamic, Dynamic, RowMajor>>())); CALL_SUBTEST_5((test_mixing_types<float>())); CALL_SUBTEST_5((test_mixed_storage()));
diff --git a/test/triangular.cpp b/test/triangular.cpp index ea4420b..55ff01e 100644 --- a/test/triangular.cpp +++ b/test/triangular.cpp
@@ -29,6 +29,14 @@ ViewType, internal::void_t<decltype(std::declval<const ViewType&>() * std::declval<typename ViewType::Scalar>())>> : std::true_type {}; +template <typename ViewType, typename = void> +struct has_structured_sum : std::false_type {}; + +template <typename ViewType> +struct has_structured_sum<ViewType, + internal::void_t<decltype(std::declval<const ViewType&>() + std::declval<const ViewType&>())>> + : std::true_type {}; + template <unsigned int Mode, typename MatrixType> void triangular_scalar_multiply(const MatrixType& m) { typedef typename MatrixType::Scalar Scalar; @@ -45,13 +53,103 @@ (triangular * s).template triangularView<Mode>().toDenseMatrix()); } +template <unsigned int Mode, typename MatrixType, bool IsSelfAdjointMode = Mode == Upper || Mode == Lower> +struct selfadjoint_structured_sum_impl { + typedef typename MatrixType::Scalar Scalar; + + static void run(const Scalar&, const MatrixType&, const MatrixType&, MatrixType&, MatrixType&) {} +}; + +template <unsigned int Mode, typename MatrixType> +struct selfadjoint_structured_sum_impl<Mode, MatrixType, true> { + typedef typename MatrixType::Scalar Scalar; + + static void run(const Scalar& s, const MatrixType& a, const MatrixType& b, MatrixType& result, + MatrixType& reference) { + if (a.rows() != a.cols()) return; + + result.setRandom(); + result.template selfadjointView<Mode>() = + s * a.template selfadjointView<Mode>() + b.template selfadjointView<Mode>(); + reference = (s * a + b).template selfadjointView<Mode>().toDenseMatrix(); + VERIFY_IS_APPROX(result, reference); + + result.setRandom(); + result.template selfadjointView<Mode>() = + a.template selfadjointView<Mode>() - s * b.template selfadjointView<Mode>(); + reference = (a - s * b).template selfadjointView<Mode>().toDenseMatrix(); + VERIFY_IS_APPROX(result, reference); + } +}; + +template <unsigned int Mode, typename MatrixType> +void triangular_structured_sum(const MatrixType& m) { + typedef typename MatrixType::Scalar Scalar; + + const Index rows = m.rows(); + const Index cols = m.cols(); + + const Scalar s = internal::random<Scalar>(); + const MatrixType a = MatrixType::Random(rows, cols); + const MatrixType b = MatrixType::Random(rows, cols); + + MatrixType result = MatrixType::Random(rows, cols); + MatrixType reference = result; + + result.template triangularView<Mode>() = s * a.template triangularView<Mode>() + b.template triangularView<Mode>(); + reference.template triangularView<Mode>() = s * a + b; + VERIFY_IS_APPROX(result, reference); + + result.setRandom(); + reference = result; + result.template triangularView<Mode>() = a.template triangularView<Mode>() - s * b.template triangularView<Mode>(); + reference.template triangularView<Mode>() = a - s * b; + VERIFY_IS_APPROX(result, reference); + + selfadjoint_structured_sum_impl<Mode, MatrixType>::run(s, a, b, result, reference); +} + +template <typename MatrixType> +void triangular_setters(const MatrixType& m) { + const Index rows = m.rows(); + const Index cols = m.cols(); + + MatrixType result = MatrixType::Random(rows, cols); + MatrixType reference = result; + + result.template triangularView<Upper>().setIdentity(); + reference.template triangularView<Upper>() = MatrixType::Identity(rows, cols); + VERIFY_IS_APPROX(result, reference); + + result.setRandom(); + reference = result; + result.template triangularView<Lower>().setIdentity(); + reference.template triangularView<Lower>() = MatrixType::Identity(rows, cols); + VERIFY_IS_APPROX(result, reference); + + result.setRandom(); + reference = result; + result.template triangularView<Upper>().setRandom(); + VERIFY_IS_APPROX(result.template triangularView<StrictlyLower>().toDenseMatrix(), + reference.template triangularView<StrictlyLower>().toDenseMatrix()); + + result.setRandom(); + reference = result; + result.template triangularView<StrictlyLower>().setRandom(); + VERIFY_IS_APPROX(result.template triangularView<Upper>().toDenseMatrix(), + reference.template triangularView<Upper>().toDenseMatrix()); +} + template <typename MatrixType> void triangular_scalar_multiply_sfinae() { + typedef decltype(std::declval<MatrixType&>().template triangularView<Upper>()) UpperView; typedef decltype(std::declval<MatrixType&>().template triangularView<Lower>()) LowerView; typedef decltype(std::declval<MatrixType&>().template triangularView<StrictlyLower>()) StrictlyLowerView; typedef decltype(std::declval<MatrixType&>().template triangularView<StrictlyUpper>()) StrictlyUpperView; typedef decltype(std::declval<MatrixType&>().template triangularView<UnitLower>()) UnitLowerView; + STATIC_CHECK((has_left_scalar_multiply<UpperView>::value)); + STATIC_CHECK((has_right_scalar_multiply<UpperView>::value)); STATIC_CHECK((has_left_scalar_multiply<LowerView>::value)); STATIC_CHECK((has_right_scalar_multiply<LowerView>::value)); STATIC_CHECK((has_left_scalar_multiply<StrictlyLowerView>::value)); @@ -60,6 +158,12 @@ STATIC_CHECK((has_right_scalar_multiply<StrictlyUpperView>::value)); STATIC_CHECK((!has_left_scalar_multiply<UnitLowerView>::value)); STATIC_CHECK((!has_right_scalar_multiply<UnitLowerView>::value)); + + STATIC_CHECK((has_structured_sum<UpperView>::value)); + STATIC_CHECK((has_structured_sum<LowerView>::value)); + STATIC_CHECK((has_structured_sum<StrictlyLowerView>::value)); + STATIC_CHECK((has_structured_sum<StrictlyUpperView>::value)); + STATIC_CHECK((!has_structured_sum<UnitLowerView>::value)); } template <typename MatrixType> @@ -92,6 +196,11 @@ typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; triangular_scalar_multiply_sfinae<MatrixType>(); + triangular_structured_sum<Upper>(m); + triangular_structured_sum<Lower>(m); + triangular_structured_sum<StrictlyUpper>(m); + triangular_structured_sum<StrictlyLower>(m); + triangular_setters(m); RealScalar largerEps = 10 * test_precision<RealScalar>(); @@ -287,6 +396,11 @@ triangular_scalar_multiply<Lower>(m1); triangular_scalar_multiply<StrictlyUpper>(m1); triangular_scalar_multiply<StrictlyLower>(m1); + triangular_structured_sum<Upper>(m1); + triangular_structured_sum<Lower>(m1); + triangular_structured_sum<StrictlyUpper>(m1); + triangular_structured_sum<StrictlyLower>(m1); + triangular_setters(m1); m1.setRandom(); m2 = m1.template triangularView<Upper>(); VERIFY(m2.isUpperTriangular());