From 7a8d0d8bc2572707c9d35006f30ea835c86954b0 Mon Sep 17 00:00:00 2001 From: sotech117 Date: Tue, 9 Apr 2024 03:14:17 -0400 Subject: first draft to generate waves --- Eigen/src/SparseCore/SparseCwiseBinaryOp.h | 722 +++++++++++++++++++++++++++++ 1 file changed, 722 insertions(+) create mode 100644 Eigen/src/SparseCore/SparseCwiseBinaryOp.h (limited to 'Eigen/src/SparseCore/SparseCwiseBinaryOp.h') diff --git a/Eigen/src/SparseCore/SparseCwiseBinaryOp.h b/Eigen/src/SparseCore/SparseCwiseBinaryOp.h new file mode 100644 index 0000000..9b0d3f9 --- /dev/null +++ b/Eigen/src/SparseCore/SparseCwiseBinaryOp.h @@ -0,0 +1,722 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2014 Gael Guennebaud +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_SPARSE_CWISE_BINARY_OP_H +#define EIGEN_SPARSE_CWISE_BINARY_OP_H + +namespace Eigen { + +// Here we have to handle 3 cases: +// 1 - sparse op dense +// 2 - dense op sparse +// 3 - sparse op sparse +// We also need to implement a 4th iterator for: +// 4 - dense op dense +// Finally, we also need to distinguish between the product and other operations : +// configuration returned mode +// 1 - sparse op dense product sparse +// generic dense +// 2 - dense op sparse product sparse +// generic dense +// 3 - sparse op sparse product sparse +// generic sparse +// 4 - dense op dense product dense +// generic dense +// +// TODO to ease compiler job, we could specialize product/quotient with a scalar +// and fallback to cwise-unary evaluator using bind1st_op and bind2nd_op. + +template +class CwiseBinaryOpImpl + : public SparseMatrixBase > +{ + public: + typedef CwiseBinaryOp Derived; + typedef SparseMatrixBase Base; + EIGEN_SPARSE_PUBLIC_INTERFACE(Derived) + CwiseBinaryOpImpl() + { + EIGEN_STATIC_ASSERT(( + (!internal::is_same::StorageKind, + typename internal::traits::StorageKind>::value) + || ((internal::evaluator::Flags&RowMajorBit) == (internal::evaluator::Flags&RowMajorBit))), + THE_STORAGE_ORDER_OF_BOTH_SIDES_MUST_MATCH); + } +}; + +namespace internal { + + +// Generic "sparse OP sparse" +template struct binary_sparse_evaluator; + +template +struct binary_evaluator, IteratorBased, IteratorBased> + : evaluator_base > +{ +protected: + typedef typename evaluator::InnerIterator LhsIterator; + typedef typename evaluator::InnerIterator RhsIterator; + typedef CwiseBinaryOp XprType; + typedef typename traits::Scalar Scalar; + typedef typename XprType::StorageIndex StorageIndex; +public: + + class InnerIterator + { + public: + + EIGEN_STRONG_INLINE InnerIterator(const binary_evaluator& aEval, Index outer) + : m_lhsIter(aEval.m_lhsImpl,outer), m_rhsIter(aEval.m_rhsImpl,outer), m_functor(aEval.m_functor) + { + this->operator++(); + } + + EIGEN_STRONG_INLINE InnerIterator& operator++() + { + if (m_lhsIter && m_rhsIter && (m_lhsIter.index() == m_rhsIter.index())) + { + m_id = m_lhsIter.index(); + m_value = m_functor(m_lhsIter.value(), m_rhsIter.value()); + ++m_lhsIter; + ++m_rhsIter; + } + else if (m_lhsIter && (!m_rhsIter || (m_lhsIter.index() < m_rhsIter.index()))) + { + m_id = m_lhsIter.index(); + m_value = m_functor(m_lhsIter.value(), Scalar(0)); + ++m_lhsIter; + } + else if (m_rhsIter && (!m_lhsIter || (m_lhsIter.index() > m_rhsIter.index()))) + { + m_id = m_rhsIter.index(); + m_value = m_functor(Scalar(0), m_rhsIter.value()); + ++m_rhsIter; + } + else + { + m_value = Scalar(0); // this is to avoid a compilation warning + m_id = -1; + } + return *this; + } + + EIGEN_STRONG_INLINE Scalar value() const { return m_value; } + + EIGEN_STRONG_INLINE StorageIndex index() const { return m_id; } + EIGEN_STRONG_INLINE Index outer() const { return m_lhsIter.outer(); } + EIGEN_STRONG_INLINE Index row() const { return Lhs::IsRowMajor ? m_lhsIter.row() : index(); } + EIGEN_STRONG_INLINE Index col() const { return Lhs::IsRowMajor ? index() : m_lhsIter.col(); } + + EIGEN_STRONG_INLINE operator bool() const { return m_id>=0; } + + protected: + LhsIterator m_lhsIter; + RhsIterator m_rhsIter; + const BinaryOp& m_functor; + Scalar m_value; + StorageIndex m_id; + }; + + + enum { + CoeffReadCost = int(evaluator::CoeffReadCost) + int(evaluator::CoeffReadCost) + int(functor_traits::Cost), + Flags = XprType::Flags + }; + + explicit binary_evaluator(const XprType& xpr) + : m_functor(xpr.functor()), + m_lhsImpl(xpr.lhs()), + m_rhsImpl(xpr.rhs()) + { + EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits::Cost); + EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); + } + + inline Index nonZerosEstimate() const { + return m_lhsImpl.nonZerosEstimate() + m_rhsImpl.nonZerosEstimate(); + } + +protected: + const BinaryOp m_functor; + evaluator m_lhsImpl; + evaluator m_rhsImpl; +}; + +// dense op sparse +template +struct binary_evaluator, IndexBased, IteratorBased> + : evaluator_base > +{ +protected: + typedef typename evaluator::InnerIterator RhsIterator; + typedef CwiseBinaryOp XprType; + typedef typename traits::Scalar Scalar; + typedef typename XprType::StorageIndex StorageIndex; +public: + + class InnerIterator + { + enum { IsRowMajor = (int(Rhs::Flags)&RowMajorBit)==RowMajorBit }; + public: + + EIGEN_STRONG_INLINE InnerIterator(const binary_evaluator& aEval, Index outer) + : m_lhsEval(aEval.m_lhsImpl), m_rhsIter(aEval.m_rhsImpl,outer), m_functor(aEval.m_functor), m_value(0), m_id(-1), m_innerSize(aEval.m_expr.rhs().innerSize()) + { + this->operator++(); + } + + EIGEN_STRONG_INLINE InnerIterator& operator++() + { + ++m_id; + if(m_id &m_lhsEval; + RhsIterator m_rhsIter; + const BinaryOp& m_functor; + Scalar m_value; + StorageIndex m_id; + StorageIndex m_innerSize; + }; + + + enum { + CoeffReadCost = int(evaluator::CoeffReadCost) + int(evaluator::CoeffReadCost) + int(functor_traits::Cost), + Flags = XprType::Flags + }; + + explicit binary_evaluator(const XprType& xpr) + : m_functor(xpr.functor()), + m_lhsImpl(xpr.lhs()), + m_rhsImpl(xpr.rhs()), + m_expr(xpr) + { + EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits::Cost); + EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); + } + + inline Index nonZerosEstimate() const { + return m_expr.size(); + } + +protected: + const BinaryOp m_functor; + evaluator m_lhsImpl; + evaluator m_rhsImpl; + const XprType &m_expr; +}; + +// sparse op dense +template +struct binary_evaluator, IteratorBased, IndexBased> + : evaluator_base > +{ +protected: + typedef typename evaluator::InnerIterator LhsIterator; + typedef CwiseBinaryOp XprType; + typedef typename traits::Scalar Scalar; + typedef typename XprType::StorageIndex StorageIndex; +public: + + class InnerIterator + { + enum { IsRowMajor = (int(Lhs::Flags)&RowMajorBit)==RowMajorBit }; + public: + + EIGEN_STRONG_INLINE InnerIterator(const binary_evaluator& aEval, Index outer) + : m_lhsIter(aEval.m_lhsImpl,outer), m_rhsEval(aEval.m_rhsImpl), m_functor(aEval.m_functor), m_value(0), m_id(-1), m_innerSize(aEval.m_expr.lhs().innerSize()) + { + this->operator++(); + } + + EIGEN_STRONG_INLINE InnerIterator& operator++() + { + ++m_id; + if(m_id &m_rhsEval; + const BinaryOp& m_functor; + Scalar m_value; + StorageIndex m_id; + StorageIndex m_innerSize; + }; + + + enum { + CoeffReadCost = int(evaluator::CoeffReadCost) + int(evaluator::CoeffReadCost) + int(functor_traits::Cost), + Flags = XprType::Flags + }; + + explicit binary_evaluator(const XprType& xpr) + : m_functor(xpr.functor()), + m_lhsImpl(xpr.lhs()), + m_rhsImpl(xpr.rhs()), + m_expr(xpr) + { + EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits::Cost); + EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); + } + + inline Index nonZerosEstimate() const { + return m_expr.size(); + } + +protected: + const BinaryOp m_functor; + evaluator m_lhsImpl; + evaluator m_rhsImpl; + const XprType &m_expr; +}; + +template::Kind, + typename RhsKind = typename evaluator_traits::Kind, + typename LhsScalar = typename traits::Scalar, + typename RhsScalar = typename traits::Scalar> struct sparse_conjunction_evaluator; + +// "sparse .* sparse" +template +struct binary_evaluator, Lhs, Rhs>, IteratorBased, IteratorBased> + : sparse_conjunction_evaluator, Lhs, Rhs> > +{ + typedef CwiseBinaryOp, Lhs, Rhs> XprType; + typedef sparse_conjunction_evaluator Base; + explicit binary_evaluator(const XprType& xpr) : Base(xpr) {} +}; +// "dense .* sparse" +template +struct binary_evaluator, Lhs, Rhs>, IndexBased, IteratorBased> + : sparse_conjunction_evaluator, Lhs, Rhs> > +{ + typedef CwiseBinaryOp, Lhs, Rhs> XprType; + typedef sparse_conjunction_evaluator Base; + explicit binary_evaluator(const XprType& xpr) : Base(xpr) {} +}; +// "sparse .* dense" +template +struct binary_evaluator, Lhs, Rhs>, IteratorBased, IndexBased> + : sparse_conjunction_evaluator, Lhs, Rhs> > +{ + typedef CwiseBinaryOp, Lhs, Rhs> XprType; + typedef sparse_conjunction_evaluator Base; + explicit binary_evaluator(const XprType& xpr) : Base(xpr) {} +}; + +// "sparse ./ dense" +template +struct binary_evaluator, Lhs, Rhs>, IteratorBased, IndexBased> + : sparse_conjunction_evaluator, Lhs, Rhs> > +{ + typedef CwiseBinaryOp, Lhs, Rhs> XprType; + typedef sparse_conjunction_evaluator Base; + explicit binary_evaluator(const XprType& xpr) : Base(xpr) {} +}; + +// "sparse && sparse" +template +struct binary_evaluator, IteratorBased, IteratorBased> + : sparse_conjunction_evaluator > +{ + typedef CwiseBinaryOp XprType; + typedef sparse_conjunction_evaluator Base; + explicit binary_evaluator(const XprType& xpr) : Base(xpr) {} +}; +// "dense && sparse" +template +struct binary_evaluator, IndexBased, IteratorBased> + : sparse_conjunction_evaluator > +{ + typedef CwiseBinaryOp XprType; + typedef sparse_conjunction_evaluator Base; + explicit binary_evaluator(const XprType& xpr) : Base(xpr) {} +}; +// "sparse && dense" +template +struct binary_evaluator, IteratorBased, IndexBased> + : sparse_conjunction_evaluator > +{ + typedef CwiseBinaryOp XprType; + typedef sparse_conjunction_evaluator Base; + explicit binary_evaluator(const XprType& xpr) : Base(xpr) {} +}; + +// "sparse ^ sparse" +template +struct sparse_conjunction_evaluator + : evaluator_base +{ +protected: + typedef typename XprType::Functor BinaryOp; + typedef typename XprType::Lhs LhsArg; + typedef typename XprType::Rhs RhsArg; + typedef typename evaluator::InnerIterator LhsIterator; + typedef typename evaluator::InnerIterator RhsIterator; + typedef typename XprType::StorageIndex StorageIndex; + typedef typename traits::Scalar Scalar; +public: + + class InnerIterator + { + public: + + EIGEN_STRONG_INLINE InnerIterator(const sparse_conjunction_evaluator& aEval, Index outer) + : m_lhsIter(aEval.m_lhsImpl,outer), m_rhsIter(aEval.m_rhsImpl,outer), m_functor(aEval.m_functor) + { + while (m_lhsIter && m_rhsIter && (m_lhsIter.index() != m_rhsIter.index())) + { + if (m_lhsIter.index() < m_rhsIter.index()) + ++m_lhsIter; + else + ++m_rhsIter; + } + } + + EIGEN_STRONG_INLINE InnerIterator& operator++() + { + ++m_lhsIter; + ++m_rhsIter; + while (m_lhsIter && m_rhsIter && (m_lhsIter.index() != m_rhsIter.index())) + { + if (m_lhsIter.index() < m_rhsIter.index()) + ++m_lhsIter; + else + ++m_rhsIter; + } + return *this; + } + + EIGEN_STRONG_INLINE Scalar value() const { return m_functor(m_lhsIter.value(), m_rhsIter.value()); } + + EIGEN_STRONG_INLINE StorageIndex index() const { return m_lhsIter.index(); } + EIGEN_STRONG_INLINE Index outer() const { return m_lhsIter.outer(); } + EIGEN_STRONG_INLINE Index row() const { return m_lhsIter.row(); } + EIGEN_STRONG_INLINE Index col() const { return m_lhsIter.col(); } + + EIGEN_STRONG_INLINE operator bool() const { return (m_lhsIter && m_rhsIter); } + + protected: + LhsIterator m_lhsIter; + RhsIterator m_rhsIter; + const BinaryOp& m_functor; + }; + + + enum { + CoeffReadCost = int(evaluator::CoeffReadCost) + int(evaluator::CoeffReadCost) + int(functor_traits::Cost), + Flags = XprType::Flags + }; + + explicit sparse_conjunction_evaluator(const XprType& xpr) + : m_functor(xpr.functor()), + m_lhsImpl(xpr.lhs()), + m_rhsImpl(xpr.rhs()) + { + EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits::Cost); + EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); + } + + inline Index nonZerosEstimate() const { + return (std::min)(m_lhsImpl.nonZerosEstimate(), m_rhsImpl.nonZerosEstimate()); + } + +protected: + const BinaryOp m_functor; + evaluator m_lhsImpl; + evaluator m_rhsImpl; +}; + +// "dense ^ sparse" +template +struct sparse_conjunction_evaluator + : evaluator_base +{ +protected: + typedef typename XprType::Functor BinaryOp; + typedef typename XprType::Lhs LhsArg; + typedef typename XprType::Rhs RhsArg; + typedef evaluator LhsEvaluator; + typedef typename evaluator::InnerIterator RhsIterator; + typedef typename XprType::StorageIndex StorageIndex; + typedef typename traits::Scalar Scalar; +public: + + class InnerIterator + { + enum { IsRowMajor = (int(RhsArg::Flags)&RowMajorBit)==RowMajorBit }; + + public: + + EIGEN_STRONG_INLINE InnerIterator(const sparse_conjunction_evaluator& aEval, Index outer) + : m_lhsEval(aEval.m_lhsImpl), m_rhsIter(aEval.m_rhsImpl,outer), m_functor(aEval.m_functor), m_outer(outer) + {} + + EIGEN_STRONG_INLINE InnerIterator& operator++() + { + ++m_rhsIter; + return *this; + } + + EIGEN_STRONG_INLINE Scalar value() const + { return m_functor(m_lhsEval.coeff(IsRowMajor?m_outer:m_rhsIter.index(),IsRowMajor?m_rhsIter.index():m_outer), m_rhsIter.value()); } + + EIGEN_STRONG_INLINE StorageIndex index() const { return m_rhsIter.index(); } + EIGEN_STRONG_INLINE Index outer() const { return m_rhsIter.outer(); } + EIGEN_STRONG_INLINE Index row() const { return m_rhsIter.row(); } + EIGEN_STRONG_INLINE Index col() const { return m_rhsIter.col(); } + + EIGEN_STRONG_INLINE operator bool() const { return m_rhsIter; } + + protected: + const LhsEvaluator &m_lhsEval; + RhsIterator m_rhsIter; + const BinaryOp& m_functor; + const Index m_outer; + }; + + + enum { + CoeffReadCost = int(evaluator::CoeffReadCost) + int(evaluator::CoeffReadCost) + int(functor_traits::Cost), + Flags = XprType::Flags + }; + + explicit sparse_conjunction_evaluator(const XprType& xpr) + : m_functor(xpr.functor()), + m_lhsImpl(xpr.lhs()), + m_rhsImpl(xpr.rhs()) + { + EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits::Cost); + EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); + } + + inline Index nonZerosEstimate() const { + return m_rhsImpl.nonZerosEstimate(); + } + +protected: + const BinaryOp m_functor; + evaluator m_lhsImpl; + evaluator m_rhsImpl; +}; + +// "sparse ^ dense" +template +struct sparse_conjunction_evaluator + : evaluator_base +{ +protected: + typedef typename XprType::Functor BinaryOp; + typedef typename XprType::Lhs LhsArg; + typedef typename XprType::Rhs RhsArg; + typedef typename evaluator::InnerIterator LhsIterator; + typedef evaluator RhsEvaluator; + typedef typename XprType::StorageIndex StorageIndex; + typedef typename traits::Scalar Scalar; +public: + + class InnerIterator + { + enum { IsRowMajor = (int(LhsArg::Flags)&RowMajorBit)==RowMajorBit }; + + public: + + EIGEN_STRONG_INLINE InnerIterator(const sparse_conjunction_evaluator& aEval, Index outer) + : m_lhsIter(aEval.m_lhsImpl,outer), m_rhsEval(aEval.m_rhsImpl), m_functor(aEval.m_functor), m_outer(outer) + {} + + EIGEN_STRONG_INLINE InnerIterator& operator++() + { + ++m_lhsIter; + return *this; + } + + EIGEN_STRONG_INLINE Scalar value() const + { return m_functor(m_lhsIter.value(), + m_rhsEval.coeff(IsRowMajor?m_outer:m_lhsIter.index(),IsRowMajor?m_lhsIter.index():m_outer)); } + + EIGEN_STRONG_INLINE StorageIndex index() const { return m_lhsIter.index(); } + EIGEN_STRONG_INLINE Index outer() const { return m_lhsIter.outer(); } + EIGEN_STRONG_INLINE Index row() const { return m_lhsIter.row(); } + EIGEN_STRONG_INLINE Index col() const { return m_lhsIter.col(); } + + EIGEN_STRONG_INLINE operator bool() const { return m_lhsIter; } + + protected: + LhsIterator m_lhsIter; + const evaluator &m_rhsEval; + const BinaryOp& m_functor; + const Index m_outer; + }; + + + enum { + CoeffReadCost = int(evaluator::CoeffReadCost) + int(evaluator::CoeffReadCost) + int(functor_traits::Cost), + Flags = XprType::Flags + }; + + explicit sparse_conjunction_evaluator(const XprType& xpr) + : m_functor(xpr.functor()), + m_lhsImpl(xpr.lhs()), + m_rhsImpl(xpr.rhs()) + { + EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits::Cost); + EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); + } + + inline Index nonZerosEstimate() const { + return m_lhsImpl.nonZerosEstimate(); + } + +protected: + const BinaryOp m_functor; + evaluator m_lhsImpl; + evaluator m_rhsImpl; +}; + +} + +/*************************************************************************** +* Implementation of SparseMatrixBase and SparseCwise functions/operators +***************************************************************************/ + +template +template +Derived& SparseMatrixBase::operator+=(const EigenBase &other) +{ + call_assignment(derived(), other.derived(), internal::add_assign_op()); + return derived(); +} + +template +template +Derived& SparseMatrixBase::operator-=(const EigenBase &other) +{ + call_assignment(derived(), other.derived(), internal::assign_op()); + return derived(); +} + +template +template +EIGEN_STRONG_INLINE Derived & +SparseMatrixBase::operator-=(const SparseMatrixBase &other) +{ + return derived() = derived() - other.derived(); +} + +template +template +EIGEN_STRONG_INLINE Derived & +SparseMatrixBase::operator+=(const SparseMatrixBase& other) +{ + return derived() = derived() + other.derived(); +} + +template +template +Derived& SparseMatrixBase::operator+=(const DiagonalBase& other) +{ + call_assignment_no_alias(derived(), other.derived(), internal::add_assign_op()); + return derived(); +} + +template +template +Derived& SparseMatrixBase::operator-=(const DiagonalBase& other) +{ + call_assignment_no_alias(derived(), other.derived(), internal::sub_assign_op()); + return derived(); +} + +template +template +EIGEN_STRONG_INLINE const typename SparseMatrixBase::template CwiseProductDenseReturnType::Type +SparseMatrixBase::cwiseProduct(const MatrixBase &other) const +{ + return typename CwiseProductDenseReturnType::Type(derived(), other.derived()); +} + +template +EIGEN_STRONG_INLINE const CwiseBinaryOp, const DenseDerived, const SparseDerived> +operator+(const MatrixBase &a, const SparseMatrixBase &b) +{ + return CwiseBinaryOp, const DenseDerived, const SparseDerived>(a.derived(), b.derived()); +} + +template +EIGEN_STRONG_INLINE const CwiseBinaryOp, const SparseDerived, const DenseDerived> +operator+(const SparseMatrixBase &a, const MatrixBase &b) +{ + return CwiseBinaryOp, const SparseDerived, const DenseDerived>(a.derived(), b.derived()); +} + +template +EIGEN_STRONG_INLINE const CwiseBinaryOp, const DenseDerived, const SparseDerived> +operator-(const MatrixBase &a, const SparseMatrixBase &b) +{ + return CwiseBinaryOp, const DenseDerived, const SparseDerived>(a.derived(), b.derived()); +} + +template +EIGEN_STRONG_INLINE const CwiseBinaryOp, const SparseDerived, const DenseDerived> +operator-(const SparseMatrixBase &a, const MatrixBase &b) +{ + return CwiseBinaryOp, const SparseDerived, const DenseDerived>(a.derived(), b.derived()); +} + +} // end namespace Eigen + +#endif // EIGEN_SPARSE_CWISE_BINARY_OP_H -- cgit v1.2.3-70-g09d2