diff --git a/Eigen/SparseCore b/Eigen/SparseCore index 69b413ec2..578469c1c 100644 --- a/Eigen/SparseCore +++ b/Eigen/SparseCore @@ -54,12 +54,12 @@ struct Sparse {}; #include "src/SparseCore/SparseProduct.h" #include "src/SparseCore/SparseDenseProduct.h" #include "src/SparseCore/SparseSelfAdjointView.h" +#include "src/SparseCore/SparseTriangularView.h" +#include "src/SparseCore/TriangularSolver.h" + #ifndef EIGEN_TEST_EVALUATORS #include "src/SparseCore/SparsePermutation.h" #include "src/SparseCore/SparseFuzzy.h" -#include "src/SparseCore/SparseTriangularView.h" - -#include "src/SparseCore/TriangularSolver.h" #endif #include "src/Core/util/ReenableStupidWarnings.h" diff --git a/Eigen/src/Core/SolveTriangular.h b/Eigen/src/Core/SolveTriangular.h index ef17f288e..0f17e3a89 100644 --- a/Eigen/src/Core/SolveTriangular.h +++ b/Eigen/src/Core/SolveTriangular.h @@ -171,10 +171,10 @@ struct triangular_solver_selector { */ template template -void TriangularView::solveInPlace(const MatrixBase& _other) const +void TriangularViewImpl::solveInPlace(const MatrixBase& _other) const { OtherDerived& other = _other.const_cast_derived(); - eigen_assert( cols() == rows() && ((Side==OnTheLeft && cols() == other.rows()) || (Side==OnTheRight && cols() == other.cols())) ); + eigen_assert( derived().cols() == derived().rows() && ((Side==OnTheLeft && derived().cols() == other.rows()) || (Side==OnTheRight && derived().cols() == other.cols())) ); eigen_assert((!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower))); enum { copy = internal::traits::Flags & RowMajorBit && OtherDerived::IsVectorAtCompileTime }; @@ -183,7 +183,7 @@ void TriangularView::solveInPlace(const MatrixBase::type, - Side, Mode>::run(nestedExpression(), otherCopy); + Side, Mode>::run(derived().nestedExpression(), otherCopy); if (copy) other = otherCopy; @@ -213,9 +213,9 @@ void TriangularView::solveInPlace(const MatrixBase template const internal::triangular_solve_retval,Other> -TriangularView::solve(const MatrixBase& other) const +TriangularViewImpl::solve(const MatrixBase& other) const { - return internal::triangular_solve_retval(*this, other.derived()); + return internal::triangular_solve_retval(derived(), other.derived()); } namespace internal { diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index 69bbaf6a9..25dab0e77 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -49,7 +49,7 @@ template class TriangularBase : public EigenBase typedef typename internal::traits::Scalar Scalar; typedef typename internal::traits::StorageKind StorageKind; typedef typename internal::traits::Index Index; - typedef typename internal::traits::DenseMatrixType DenseMatrixType; + typedef typename internal::traits::FullMatrixType DenseMatrixType; typedef DenseMatrixType DenseType; typedef Derived const& Nested; @@ -172,8 +172,8 @@ struct traits > : traits typedef typename nested::type MatrixTypeNested; typedef typename remove_reference::type MatrixTypeNestedNonRef; typedef typename remove_all::type MatrixTypeNestedCleaned; + typedef typename MatrixType::PlainObject FullMatrixType; typedef MatrixType ExpressionType; - typedef typename MatrixType::PlainObject DenseMatrixType; enum { Mode = _Mode, Flags = (MatrixTypeNestedCleaned::Flags & (HereditaryBits | LvalueBit) & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit))) | Mode @@ -185,22 +185,23 @@ struct traits > : traits }; } +#ifndef EIGEN_TEST_EVALUATORS template struct TriangularProduct; +#endif + +template class TriangularViewImpl; template class TriangularView - : public TriangularBase > + : public TriangularViewImpl<_MatrixType, _Mode, typename internal::traits<_MatrixType>::StorageKind > { public: - typedef TriangularBase Base; + typedef TriangularViewImpl<_MatrixType, _Mode, typename internal::traits<_MatrixType>::StorageKind > Base; typedef typename internal::traits::Scalar Scalar; - typedef _MatrixType MatrixType; - typedef typename internal::traits::DenseMatrixType DenseMatrixType; - typedef DenseMatrixType PlainObject; protected: typedef typename internal::traits::MatrixTypeNested MatrixTypeNested; @@ -210,8 +211,6 @@ template class TriangularView typedef typename internal::remove_all::type MatrixConjugateReturnType; public: - using Base::evalToLazy; - typedef typename internal::traits::StorageKind StorageKind; typedef typename internal::traits::Index Index; @@ -228,110 +227,19 @@ template class TriangularView EIGEN_DEVICE_FUNC inline TriangularView(const MatrixType& matrix) : m_matrix(matrix) {} + + using Base::operator=; EIGEN_DEVICE_FUNC inline Index rows() const { return m_matrix.rows(); } EIGEN_DEVICE_FUNC inline Index cols() const { return m_matrix.cols(); } - EIGEN_DEVICE_FUNC - inline Index outerStride() const { return m_matrix.outerStride(); } - EIGEN_DEVICE_FUNC - inline Index innerStride() const { return m_matrix.innerStride(); } - -#ifdef EIGEN_TEST_EVALUATORS - - /** \sa MatrixBase::operator+=() */ - template - EIGEN_DEVICE_FUNC - TriangularView& operator+=(const DenseBase& other) { - internal::call_assignment_no_alias(*this, other.derived(), internal::add_assign_op()); - return *this; - } - /** \sa MatrixBase::operator-=() */ - template - EIGEN_DEVICE_FUNC - TriangularView& operator-=(const DenseBase& other) { - internal::call_assignment_no_alias(*this, other.derived(), internal::sub_assign_op()); - return *this; - } - -#else - /** \sa MatrixBase::operator+=() */ - template - EIGEN_DEVICE_FUNC - TriangularView& operator+=(const DenseBase& other) { return *this = m_matrix + other.derived(); } - /** \sa MatrixBase::operator-=() */ - template - EIGEN_DEVICE_FUNC - TriangularView& operator-=(const DenseBase& other) { return *this = m_matrix - other.derived(); } -#endif - /** \sa MatrixBase::operator*=() */ - EIGEN_DEVICE_FUNC - TriangularView& operator*=(const typename internal::traits::Scalar& other) { return *this = m_matrix * other; } - /** \sa MatrixBase::operator/=() */ - EIGEN_DEVICE_FUNC - TriangularView& operator/=(const typename internal::traits::Scalar& other) { return *this = m_matrix / other; } - - /** \sa MatrixBase::fill() */ - EIGEN_DEVICE_FUNC - void fill(const Scalar& value) { setConstant(value); } - /** \sa MatrixBase::setConstant() */ - EIGEN_DEVICE_FUNC - TriangularView& setConstant(const Scalar& value) - { return *this = MatrixType::Constant(rows(), cols(), value); } - /** \sa MatrixBase::setZero() */ - EIGEN_DEVICE_FUNC - TriangularView& setZero() { return setConstant(Scalar(0)); } - /** \sa MatrixBase::setOnes() */ - EIGEN_DEVICE_FUNC - TriangularView& 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 m_matrix.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) - { - Base::check_coordinates_internal(row, col); - return m_matrix.const_cast_derived().coeffRef(row, col); - } EIGEN_DEVICE_FUNC const MatrixTypeNestedCleaned& nestedExpression() const { return m_matrix; } EIGEN_DEVICE_FUNC MatrixTypeNestedCleaned& nestedExpression() { return *const_cast(&m_matrix); } - /** Assigns a triangular matrix to a triangular part of a dense matrix */ - template - EIGEN_DEVICE_FUNC - TriangularView& operator=(const TriangularBase& other); - - template - EIGEN_DEVICE_FUNC - TriangularView& operator=(const MatrixBase& other); - - EIGEN_DEVICE_FUNC - TriangularView& operator=(const TriangularView& other) - { return *this = other.nestedExpression(); } - - template - EIGEN_DEVICE_FUNC - void lazyAssign(const TriangularBase& other); - - template - EIGEN_DEVICE_FUNC - void lazyAssign(const MatrixBase& other); - /** \sa MatrixBase::conjugate() */ EIGEN_DEVICE_FUNC inline TriangularView conjugate() @@ -360,78 +268,16 @@ template class TriangularView return m_matrix.transpose(); } -#ifdef EIGEN_TEST_EVALUATORS - - /** Efficient triangular matrix times vector/matrix product */ - template - EIGEN_DEVICE_FUNC - const Product - operator*(const MatrixBase& rhs) const - { - return Product(*this, rhs.derived()); - } - - /** Efficient vector/matrix times triangular matrix product */ - template friend - EIGEN_DEVICE_FUNC - const Product - operator*(const MatrixBase& lhs, const TriangularView& rhs) - { - return Product(lhs.derived(),rhs); - } - -#else // EIGEN_TEST_EVALUATORS - /** Efficient triangular matrix times vector/matrix product */ - template - EIGEN_DEVICE_FUNC - TriangularProduct - operator*(const MatrixBase& rhs) const - { - return TriangularProduct - - (m_matrix, rhs.derived()); - } - - /** Efficient vector/matrix times triangular matrix product */ - template friend - EIGEN_DEVICE_FUNC - TriangularProduct - operator*(const MatrixBase& lhs, const TriangularView& rhs) - { - return TriangularProduct - - (lhs.derived(),rhs.m_matrix); - } -#endif - - template - EIGEN_DEVICE_FUNC - inline const internal::triangular_solve_retval - solve(const MatrixBase& other) const; - - template - EIGEN_DEVICE_FUNC - void solveInPlace(const MatrixBase& other) const; - #ifdef EIGEN_TEST_EVALUATORS template EIGEN_DEVICE_FUNC inline const Solve solve(const MatrixBase& other) const { return Solve(*this, other.derived()); } -#else // EIGEN_TEST_EVALUATORS - template - EIGEN_DEVICE_FUNC - inline const internal::triangular_solve_retval - solve(const MatrixBase& other) const - { return solve(other); } + + using Base::solve; #endif // EIGEN_TEST_EVALUATORS - template - EIGEN_DEVICE_FUNC - void solveInPlace(const MatrixBase& other) const - { return solveInPlace(other); } - EIGEN_DEVICE_FUNC const SelfAdjointView selfadjointView() const { @@ -445,30 +291,6 @@ template class TriangularView return SelfAdjointView(m_matrix); } - template - EIGEN_DEVICE_FUNC - void swap(TriangularBase const & other) - { - #ifdef EIGEN_TEST_EVALUATORS - call_assignment(*this, other.const_cast_derived(), internal::swap_assign_op()); - #else - TriangularView,Mode>(const_cast(m_matrix)).lazyAssign(other.const_cast_derived().nestedExpression()); - #endif - } - - // TODO: this overload is ambiguous and it should be deprecated (Gael) - template - EIGEN_DEVICE_FUNC - void swap(MatrixBase const & other) - { - #ifdef EIGEN_TEST_EVALUATORS - call_assignment(*this, other.const_cast_derived(), internal::swap_assign_op()); - #else - SwapWrapper swaper(const_cast(m_matrix)); - TriangularView,Mode>(swaper).lazyAssign(other.derived()); - #endif - } - EIGEN_DEVICE_FUNC Scalar determinant() const { @@ -479,13 +301,227 @@ template class TriangularView else return m_matrix.diagonal().prod(); } + + protected: + + MatrixTypeNested m_matrix; +}; + +template class TriangularViewImpl<_MatrixType,_Mode,Dense> + : public TriangularBase > +{ + public: + + typedef TriangularView<_MatrixType, _Mode> TriangularViewType; + typedef TriangularBase Base; + typedef typename internal::traits::Scalar Scalar; + + typedef _MatrixType MatrixType; + typedef typename MatrixType::PlainObject DenseMatrixType; + typedef DenseMatrixType PlainObject; + + public: + using Base::evalToLazy; + using Base::derived; + + typedef typename internal::traits::StorageKind StorageKind; + typedef typename internal::traits::Index Index; + + enum { + Mode = _Mode, + Flags = internal::traits::Flags + }; + + EIGEN_DEVICE_FUNC + inline Index outerStride() const { return derived().nestedExpression().outerStride(); } + EIGEN_DEVICE_FUNC + inline Index innerStride() const { return derived().nestedExpression().innerStride(); } + +#ifdef EIGEN_TEST_EVALUATORS + + /** \sa MatrixBase::operator+=() */ + template + EIGEN_DEVICE_FUNC + TriangularViewType& operator+=(const DenseBase& other) { + internal::call_assignment_no_alias(derived(), other.derived(), internal::add_assign_op()); + return derived(); + } + /** \sa MatrixBase::operator-=() */ + template + EIGEN_DEVICE_FUNC + TriangularViewType& operator-=(const DenseBase& other) { + internal::call_assignment_no_alias(derived(), other.derived(), internal::sub_assign_op()); + return derived(); + } + +#else + /** \sa MatrixBase::operator+=() */ + template + EIGEN_DEVICE_FUNC + TriangularViewType& operator+=(const DenseBase& other) { return *this = derived().nestedExpression() + other.derived(); } + /** \sa MatrixBase::operator-=() */ + template + EIGEN_DEVICE_FUNC + TriangularViewType& operator-=(const DenseBase& other) { return *this = derived().nestedExpression() - other.derived(); } +#endif + /** \sa MatrixBase::operator*=() */ + EIGEN_DEVICE_FUNC + TriangularViewType& operator*=(const typename internal::traits::Scalar& other) { return *this = derived().nestedExpression() * other; } + /** \sa MatrixBase::operator/=() */ + EIGEN_DEVICE_FUNC + TriangularViewType& operator/=(const typename internal::traits::Scalar& other) { 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) + { + Base::check_coordinates_internal(row, col); + return derived().nestedExpression().const_cast_derived().coeffRef(row, col); + } + + /** Assigns a triangular matrix to a triangular part of a dense matrix */ + template + EIGEN_DEVICE_FUNC + TriangularViewType& operator=(const TriangularBase& other); + + template + EIGEN_DEVICE_FUNC + TriangularViewType& operator=(const MatrixBase& other); + + EIGEN_DEVICE_FUNC + TriangularViewType& operator=(const TriangularViewImpl& other) + { return *this = other.nestedExpression(); } + + template + EIGEN_DEVICE_FUNC + void lazyAssign(const TriangularBase& other); + + template + EIGEN_DEVICE_FUNC + void lazyAssign(const MatrixBase& other); + +#ifdef EIGEN_TEST_EVALUATORS + + /** Efficient triangular matrix times vector/matrix product */ + template + EIGEN_DEVICE_FUNC + const Product + operator*(const MatrixBase& rhs) const + { + return Product(derived(), rhs.derived()); + } + + /** Efficient vector/matrix times triangular matrix product */ + template friend + EIGEN_DEVICE_FUNC + const Product + operator*(const MatrixBase& lhs, const TriangularViewImpl& rhs) + { + return Product(lhs.derived(),rhs.derived()); + } + +#else // EIGEN_TEST_EVALUATORS + /** Efficient triangular matrix times vector/matrix product */ + template + EIGEN_DEVICE_FUNC + TriangularProduct + operator*(const MatrixBase& rhs) const + { + return TriangularProduct + + (derived().nestedExpression(), rhs.derived()); + } + + /** Efficient vector/matrix times triangular matrix product */ + template friend + EIGEN_DEVICE_FUNC + TriangularProduct + operator*(const MatrixBase& lhs, const TriangularVieImplw& rhs) + { + return TriangularProduct + + (lhs.derived(),rhs.nestedExpression()); + } +#endif + + template + EIGEN_DEVICE_FUNC + inline const internal::triangular_solve_retval + solve(const MatrixBase& other) const; + + template + EIGEN_DEVICE_FUNC + void solveInPlace(const MatrixBase& other) const; + +#ifndef EIGEN_TEST_EVALUATORS + template + EIGEN_DEVICE_FUNC + inline const internal::triangular_solve_retval + solve(const MatrixBase& other) const + { return solve(other); } +#endif // EIGEN_TEST_EVALUATORS + + template + EIGEN_DEVICE_FUNC + void solveInPlace(const MatrixBase& other) const + { return solveInPlace(other); } + + template + EIGEN_DEVICE_FUNC + void swap(TriangularBase const & other) + { + #ifdef EIGEN_TEST_EVALUATORS + call_assignment(derived(), other.const_cast_derived(), internal::swap_assign_op()); + #else + TriangularView,Mode>(const_cast(derived().nestedExpression())).lazyAssign(other.const_cast_derived().nestedExpression()); + #endif + } + + // TODO: this overload is ambiguous and it should be deprecated (Gael) + template + EIGEN_DEVICE_FUNC + void swap(MatrixBase const & other) + { + #ifdef EIGEN_TEST_EVALUATORS + call_assignment(derived(), other.const_cast_derived(), internal::swap_assign_op()); + #else + SwapWrapper swaper(const_cast(derived().nestedExpression())); + TriangularView,Mode>(swaper).lazyAssign(other.derived()); + #endif + } #ifndef EIGEN_TEST_EVALUATORS // TODO simplify the following: template EIGEN_DEVICE_FUNC - EIGEN_STRONG_INLINE TriangularView& operator=(const ProductBase& other) + EIGEN_STRONG_INLINE TriangularViewType& operator=(const ProductBase& other) { setZero(); return assignProduct(other,1); @@ -493,14 +529,14 @@ template class TriangularView template EIGEN_DEVICE_FUNC - EIGEN_STRONG_INLINE TriangularView& operator+=(const ProductBase& other) + EIGEN_STRONG_INLINE TriangularViewType& operator+=(const ProductBase& other) { return assignProduct(other,1); } template EIGEN_DEVICE_FUNC - EIGEN_STRONG_INLINE TriangularView& operator-=(const ProductBase& other) + EIGEN_STRONG_INLINE TriangularViewType& operator-=(const ProductBase& other) { return assignProduct(other,-1); } @@ -508,7 +544,7 @@ template class TriangularView template EIGEN_DEVICE_FUNC - EIGEN_STRONG_INLINE TriangularView& operator=(const ScaledProduct& other) + EIGEN_STRONG_INLINE TriangularViewType& operator=(const ScaledProduct& other) { setZero(); return assignProduct(other,other.alpha()); @@ -516,14 +552,14 @@ template class TriangularView template EIGEN_DEVICE_FUNC - EIGEN_STRONG_INLINE TriangularView& operator+=(const ScaledProduct& other) + EIGEN_STRONG_INLINE TriangularViewType& operator+=(const ScaledProduct& other) { return assignProduct(other,other.alpha()); } template EIGEN_DEVICE_FUNC - EIGEN_STRONG_INLINE TriangularView& operator-=(const ScaledProduct& other) + EIGEN_STRONG_INLINE TriangularViewType& operator-=(const ScaledProduct& other) { return assignProduct(other,-other.alpha()); } @@ -542,16 +578,15 @@ template class TriangularView template EIGEN_DEVICE_FUNC - EIGEN_STRONG_INLINE TriangularView& _assignProduct(const ProductType& prod, const Scalar& alpha); + EIGEN_STRONG_INLINE TriangularViewType& _assignProduct(const ProductType& prod, const Scalar& alpha); protected: #else protected: template EIGEN_DEVICE_FUNC - EIGEN_STRONG_INLINE TriangularView& assignProduct(const ProductBase& prod, const Scalar& alpha); + EIGEN_STRONG_INLINE TriangularViewType& assignProduct(const ProductBase& prod, const Scalar& alpha); #endif - MatrixTypeNested m_matrix; }; /*************************************************************************** @@ -736,18 +771,18 @@ struct triangular_assignment_selector template inline TriangularView& -TriangularView::operator=(const MatrixBase& other) +TriangularViewImpl::operator=(const MatrixBase& other) { - internal::call_assignment_no_alias(*this, other.derived(), internal::assign_op()); - return *this; + internal::call_assignment_no_alias(derived(), other.derived(), internal::assign_op()); + return derived(); } // FIXME should we keep that possibility template template -void TriangularView::lazyAssign(const MatrixBase& other) +void TriangularViewImpl::lazyAssign(const MatrixBase& other) { - internal::call_assignment(this->noalias(), other.template triangularView()); + internal::call_assignment(derived().noalias(), other.template triangularView()); } @@ -755,19 +790,19 @@ void TriangularView::lazyAssign(const MatrixBase template template inline TriangularView& -TriangularView::operator=(const TriangularBase& other) +TriangularViewImpl::operator=(const TriangularBase& other) { eigen_assert(Mode == int(OtherDerived::Mode)); - internal::call_assignment(*this, other.derived()); - return *this; + internal::call_assignment(derived(), other.derived()); + return derived(); } template template -void TriangularView::lazyAssign(const TriangularBase& other) +void TriangularViewImpl::lazyAssign(const TriangularBase& other) { eigen_assert(Mode == int(OtherDerived::Mode)); - internal::call_assignment(this->noalias(), other.derived()); + internal::call_assignment(derived().noalias(), other.derived()); } #else @@ -776,7 +811,7 @@ void TriangularView::lazyAssign(const TriangularBase template inline TriangularView& -TriangularView::operator=(const MatrixBase& other) +TriangularViewImpl::operator=(const MatrixBase& other) { if(OtherDerived::Flags & EvalBeforeAssigningBit) { @@ -786,13 +821,13 @@ TriangularView::operator=(const MatrixBase& othe } else lazyAssign(other.derived()); - return *this; + return derived(); } // FIXME should we keep that possibility template template -void TriangularView::lazyAssign(const MatrixBase& other) +void TriangularViewImp::lazyAssign(const MatrixBase& other) { enum { unroll = MatrixType::SizeAtCompileTime != Dynamic @@ -813,7 +848,7 @@ void TriangularView::lazyAssign(const MatrixBase template template inline TriangularView& -TriangularView::operator=(const TriangularBase& other) +TriangularViewImpl::operator=(const TriangularBase& other) { eigen_assert(Mode == int(OtherDerived::Mode)); if(internal::traits::Flags & EvalBeforeAssigningBit) @@ -824,12 +859,12 @@ TriangularView::operator=(const TriangularBase& } else lazyAssign(other.derived().nestedExpression()); - return *this; + return derived(); } template template -void TriangularView::lazyAssign(const TriangularBase& other) +void TriangularViewImpl::lazyAssign(const TriangularBase& other) { enum { unroll = MatrixType::SizeAtCompileTime != Dynamic @@ -1000,7 +1035,7 @@ template struct evaluator_traits > { typedef typename storage_kind_to_evaluator_kind::Kind Kind; - typedef TriangularShape Shape; + typedef typename glue_shapes::Shape, TriangularShape>::type Shape; // 1 if assignment A = B assumes aliasing when B is of type T and thus B needs to be evaluated into a // temporary; 0 if not. diff --git a/Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h b/Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h index b28f07bdc..a9d8352a9 100644 --- a/Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h +++ b/Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h @@ -264,22 +264,23 @@ struct general_product_to_triangular_selector #ifdef EIGEN_TEST_EVALUATORS template template -TriangularView& TriangularView::_assignProduct(const ProductType& prod, const Scalar& alpha) +TriangularView& TriangularViewImpl::_assignProduct(const ProductType& prod, const Scalar& alpha) { - eigen_assert(m_matrix.rows() == prod.rows() && m_matrix.cols() == prod.cols()); + eigen_assert(derived().nestedExpression().rows() == prod.rows() && derived().cols() == prod.cols()); - general_product_to_triangular_selector::InnerSize==1>::run(m_matrix.const_cast_derived(), prod, alpha); + general_product_to_triangular_selector::InnerSize==1>::run(derived().nestedExpression().const_cast_derived(), prod, alpha); return *this; } #else template template -TriangularView& TriangularView::assignProduct(const ProductBase& prod, const Scalar& alpha) +TriangularView& TriangularViewImpl::assignProduct(const ProductBase& prod, const Scalar& alpha) { - eigen_assert(m_matrix.rows() == prod.rows() && m_matrix.cols() == prod.cols()); + eigen_assert(derived().rows() == prod.rows() && derived().cols() == prod.cols()); - general_product_to_triangular_selector::run(m_matrix.const_cast_derived(), prod.derived(), alpha); + general_product_to_triangular_selector + ::run(derived().nestedExpression().const_cast_derived(), prod.derived(), alpha); return *this; } diff --git a/Eigen/src/Core/products/TriangularMatrixMatrix.h b/Eigen/src/Core/products/TriangularMatrixMatrix.h index d93277cb8..fda6e2486 100644 --- a/Eigen/src/Core/products/TriangularMatrixMatrix.h +++ b/Eigen/src/Core/products/TriangularMatrixMatrix.h @@ -368,10 +368,12 @@ EIGEN_DONT_INLINE void product_triangular_matrix_matrix struct traits > : traits, Lhs, Rhs> > {}; +#endif } // end namespace internal diff --git a/Eigen/src/Core/products/TriangularMatrixVector.h b/Eigen/src/Core/products/TriangularMatrixVector.h index 399de3092..19167c232 100644 --- a/Eigen/src/Core/products/TriangularMatrixVector.h +++ b/Eigen/src/Core/products/TriangularMatrixVector.h @@ -157,6 +157,8 @@ EIGEN_DONT_INLINE void triangular_matrix_vector_product struct traits > : traits, Lhs, Rhs> > @@ -166,7 +168,7 @@ template struct traits > : traits, Lhs, Rhs> > {}; - +#endif template struct trmv_selector; diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 4b55b3a9a..3ab8d0ed3 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -450,13 +450,13 @@ struct MatrixXpr {}; struct ArrayXpr {}; // An evaluator must define its shape. By default, it can be one of the following: -struct DenseShape { static std::string debugName() { return "DenseShape"; } }; -struct DiagonalShape { static std::string debugName() { return "DiagonalShape"; } }; -struct BandShape { static std::string debugName() { return "BandShape"; } }; -struct TriangularShape { static std::string debugName() { return "TriangularShape"; } }; -struct SelfAdjointShape { static std::string debugName() { return "SelfAdjointShape"; } }; -struct PermutationShape { static std::string debugName() { return "PermutationShape"; } }; -struct SparseShape { static std::string debugName() { return "SparseShape"; } }; +struct DenseShape { static std::string debugName() { return "DenseShape"; } }; +struct DiagonalShape { static std::string debugName() { return "DiagonalShape"; } }; +struct BandShape { static std::string debugName() { return "BandShape"; } }; +struct TriangularShape { static std::string debugName() { return "TriangularShape"; } }; +struct SelfAdjointShape { static std::string debugName() { return "SelfAdjointShape"; } }; +struct PermutationShape { static std::string debugName() { return "PermutationShape"; } }; +struct SparseShape { static std::string debugName() { return "SparseShape"; } }; } // end namespace Eigen diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h index 70f2c566f..7779fb3fa 100644 --- a/Eigen/src/Core/util/XprHelper.h +++ b/Eigen/src/Core/util/XprHelper.h @@ -638,6 +638,9 @@ template struct is_diagonal > template struct is_diagonal > { enum { ret = true }; }; +template struct glue_shapes; +template<> struct glue_shapes { typedef TriangularShape type; }; + } // end namespace internal // we require Lhs and Rhs to have the same scalar type. Currently there is no example of a binary functor diff --git a/Eigen/src/SparseCholesky/SimplicialCholesky.h b/Eigen/src/SparseCholesky/SimplicialCholesky.h index f41d7e010..21c7e0b39 100644 --- a/Eigen/src/SparseCholesky/SimplicialCholesky.h +++ b/Eigen/src/SparseCholesky/SimplicialCholesky.h @@ -253,8 +253,8 @@ template struct traits CholMatrixType; - typedef SparseTriangularView MatrixL; - typedef SparseTriangularView MatrixU; + typedef TriangularView MatrixL; + typedef TriangularView MatrixU; static inline MatrixL getL(const MatrixType& m) { return m; } static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); } }; @@ -266,8 +266,8 @@ template struct traits CholMatrixType; - typedef SparseTriangularView MatrixL; - typedef SparseTriangularView MatrixU; + typedef TriangularView MatrixL; + typedef TriangularView MatrixU; static inline MatrixL getL(const MatrixType& m) { return m; } static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); } }; diff --git a/Eigen/src/SparseCore/MappedSparseMatrix.h b/Eigen/src/SparseCore/MappedSparseMatrix.h index ab1a266a9..9205b906f 100644 --- a/Eigen/src/SparseCore/MappedSparseMatrix.h +++ b/Eigen/src/SparseCore/MappedSparseMatrix.h @@ -176,6 +176,34 @@ class MappedSparseMatrix::ReverseInnerIterator const Index m_end; }; +#ifdef EIGEN_ENABLE_EVALUATORS +namespace internal { + +template +struct evaluator > + : evaluator_base > +{ + typedef MappedSparseMatrix<_Scalar,_Options,_Index> MappedSparseMatrixType; + typedef typename MappedSparseMatrixType::InnerIterator InnerIterator; + typedef typename MappedSparseMatrixType::ReverseInnerIterator ReverseInnerIterator; + + enum { + CoeffReadCost = NumTraits<_Scalar>::ReadCost, + Flags = MappedSparseMatrixType::Flags + }; + + evaluator() : m_matrix(0) {} + evaluator(const MappedSparseMatrixType &mat) : m_matrix(&mat) {} + + operator MappedSparseMatrixType&() { return m_matrix->const_cast_derived(); } + operator const MappedSparseMatrixType&() const { return *m_matrix; } + + const MappedSparseMatrixType *m_matrix; +}; + +} +#endif + } // end namespace Eigen #endif // EIGEN_MAPPED_SPARSEMATRIX_H diff --git a/Eigen/src/SparseCore/SparseMatrixBase.h b/Eigen/src/SparseCore/SparseMatrixBase.h index 361a66067..4c5563755 100644 --- a/Eigen/src/SparseCore/SparseMatrixBase.h +++ b/Eigen/src/SparseCore/SparseMatrixBase.h @@ -333,7 +333,7 @@ template class SparseMatrixBase : public EigenBase Derived& operator*=(const SparseMatrixBase& other); template - inline const SparseTriangularView triangularView() const; + inline const TriangularView triangularView() const; template inline const SparseSelfAdjointView selfadjointView() const; template inline SparseSelfAdjointView selfadjointView(); diff --git a/Eigen/src/SparseCore/SparseSelfAdjointView.h b/Eigen/src/SparseCore/SparseSelfAdjointView.h index 8bd836883..ff51fc435 100644 --- a/Eigen/src/SparseCore/SparseSelfAdjointView.h +++ b/Eigen/src/SparseCore/SparseSelfAdjointView.h @@ -392,8 +392,6 @@ inline void sparse_selfadjoint_time_dense_product(const SparseLhsType& lhs, cons res.row(j) += i.value() * rhs.row(j); } } - -struct SparseSelfAdjointShape { static std::string debugName() { return "SparseSelfAdjointShape"; } }; // TODO currently a selfadjoint expression has the form SelfAdjointView<.,.> // in the future selfadjoint-ness should be defined by the expression traits diff --git a/Eigen/src/SparseCore/SparseTriangularView.h b/Eigen/src/SparseCore/SparseTriangularView.h index 333127b78..ad184208b 100644 --- a/Eigen/src/SparseCore/SparseTriangularView.h +++ b/Eigen/src/SparseCore/SparseTriangularView.h @@ -15,15 +15,15 @@ namespace Eigen { namespace internal { -template -struct traits > -: public traits -{}; +// template +// struct traits > +// : public traits +// {}; } // namespace internal -template class SparseTriangularView - : public SparseMatrixBase > +template class TriangularViewImpl + : public SparseMatrixBase > { enum { SkipFirst = ((Mode&Lower) && !(MatrixType::Flags&RowMajorBit)) || ((Mode&Upper) && (MatrixType::Flags&RowMajorBit)), @@ -31,45 +31,51 @@ template class SparseTriangularView SkipDiag = (Mode&ZeroDiag) ? 1 : 0, HasUnitDiag = (Mode&UnitDiag) ? 1 : 0 }; + + typedef TriangularView TriangularViewType; + +protected: + // dummy solve function to make TriangularView happy. + void solve() const; public: - EIGEN_SPARSE_PUBLIC_INTERFACE(SparseTriangularView) - + EIGEN_SPARSE_PUBLIC_INTERFACE(TriangularViewType) + class InnerIterator; class ReverseInnerIterator; - inline Index rows() const { return m_matrix.rows(); } - inline Index cols() const { return m_matrix.cols(); } - typedef typename MatrixType::Nested MatrixTypeNested; typedef typename internal::remove_reference::type MatrixTypeNestedNonRef; typedef typename internal::remove_all::type MatrixTypeNestedCleaned; - inline SparseTriangularView(const MatrixType& matrix) : m_matrix(matrix) {} - - /** \internal */ - inline const MatrixTypeNestedCleaned& nestedExpression() const { return m_matrix; } - +#ifndef EIGEN_TEST_EVALUATORS template typename internal::plain_matrix_type_column_major::type solve(const MatrixBase& other) const; +#else // EIGEN_TEST_EVALUATORS + template + EIGEN_DEVICE_FUNC + EIGEN_STRONG_INLINE void _solve_impl(const RhsType &rhs, DstType &dst) const { + if(!(internal::is_same::value && internal::extract_data(dst) == internal::extract_data(rhs))) + dst = rhs; + this->template solveInPlace(dst); + } +#endif // EIGEN_TEST_EVALUATORS template void solveInPlace(MatrixBase& other) const; template void solveInPlace(SparseMatrixBase& other) const; - - protected: - MatrixTypeNested m_matrix; + }; -template -class SparseTriangularView::InnerIterator : public MatrixTypeNestedCleaned::InnerIterator +template +class TriangularViewImpl::InnerIterator : public MatrixTypeNestedCleaned::InnerIterator { typedef typename MatrixTypeNestedCleaned::InnerIterator Base; - typedef typename SparseTriangularView::Index Index; + typedef typename TriangularViewType::Index Index; public: - EIGEN_STRONG_INLINE InnerIterator(const SparseTriangularView& view, Index outer) + EIGEN_STRONG_INLINE InnerIterator(const TriangularViewImpl& view, Index outer) : Base(view.nestedExpression(), outer), m_returnOne(false) { if(SkipFirst) @@ -132,14 +138,14 @@ class SparseTriangularView::InnerIterator : public MatrixTypeNe bool m_returnOne; }; -template -class SparseTriangularView::ReverseInnerIterator : public MatrixTypeNestedCleaned::ReverseInnerIterator +template +class TriangularViewImpl::ReverseInnerIterator : public MatrixTypeNestedCleaned::ReverseInnerIterator { typedef typename MatrixTypeNestedCleaned::ReverseInnerIterator Base; - typedef typename SparseTriangularView::Index Index; + typedef typename TriangularViewImpl::Index Index; public: - EIGEN_STRONG_INLINE ReverseInnerIterator(const SparseTriangularView& view, Index outer) + EIGEN_STRONG_INLINE ReverseInnerIterator(const TriangularViewType& view, Index outer) : Base(view.nestedExpression(), outer) { eigen_assert((!HasUnitDiag) && "ReverseInnerIterator does not support yet triangular views with a unit diagonal"); @@ -168,7 +174,7 @@ class SparseTriangularView::ReverseInnerIterator : public Matri template template -inline const SparseTriangularView +inline const TriangularView SparseMatrixBase::triangularView() const { return derived(); diff --git a/Eigen/src/SparseCore/SparseUtil.h b/Eigen/src/SparseCore/SparseUtil.h index 6e1db0ce8..686be6c92 100644 --- a/Eigen/src/SparseCore/SparseUtil.h +++ b/Eigen/src/SparseCore/SparseUtil.h @@ -94,7 +94,6 @@ template class Dynamic template class SparseVector; template class MappedSparseMatrix; -template class SparseTriangularView; template class SparseSelfAdjointView; template class SparseDiagonalProduct; template class SparseView; @@ -163,6 +162,12 @@ struct generic_xpr_base typedef SparseMatrixBase type; }; +struct SparseTriangularShape { static std::string debugName() { return "SparseTriangularShape"; } }; +struct SparseSelfAdjointShape { static std::string debugName() { return "SparseSelfAdjointShape"; } }; + +template<> struct glue_shapes { typedef SparseSelfAdjointShape type; }; +template<> struct glue_shapes { typedef SparseTriangularShape type; }; + } // end namespace internal /** \ingroup SparseCore_Module diff --git a/Eigen/src/SparseCore/TriangularSolver.h b/Eigen/src/SparseCore/TriangularSolver.h index dd55522a7..012a1bb75 100644 --- a/Eigen/src/SparseCore/TriangularSolver.h +++ b/Eigen/src/SparseCore/TriangularSolver.h @@ -23,6 +23,7 @@ template::Flags) & RowMajorBit> struct sparse_solve_triangular_selector; +#ifndef EIGEN_TEST_EVALUATORS // forward substitution, row-major template struct sparse_solve_triangular_selector @@ -163,13 +164,166 @@ struct sparse_solve_triangular_selector } }; +#else // EIGEN_TEST_EVALUATORS + +// forward substitution, row-major +template +struct sparse_solve_triangular_selector +{ + typedef typename Rhs::Scalar Scalar; + typedef typename Lhs::Index Index; + typedef typename evaluator::type LhsEval; + typedef typename evaluator::InnerIterator LhsIterator; + static void run(const Lhs& lhs, Rhs& other) + { + LhsEval lhsEval(lhs); + for(Index col=0 ; col +struct sparse_solve_triangular_selector +{ + typedef typename Rhs::Scalar Scalar; + typedef typename Lhs::Index Index; + typedef typename evaluator::type LhsEval; + typedef typename evaluator::InnerIterator LhsIterator; + static void run(const Lhs& lhs, Rhs& other) + { + LhsEval lhsEval(lhs); + for(Index col=0 ; col=0 ; --i) + { + Scalar tmp = other.coeff(i,col); + Scalar l_ii = 0; + LhsIterator it(lhsEval, i); + while(it && it.index() +struct sparse_solve_triangular_selector +{ + typedef typename Rhs::Scalar Scalar; + typedef typename Lhs::Index Index; + typedef typename evaluator::type LhsEval; + typedef typename evaluator::InnerIterator LhsIterator; + static void run(const Lhs& lhs, Rhs& other) + { + LhsEval lhsEval(lhs); + for(Index col=0 ; col +struct sparse_solve_triangular_selector +{ + typedef typename Rhs::Scalar Scalar; + typedef typename Lhs::Index Index; + typedef typename evaluator::type LhsEval; + typedef typename evaluator::InnerIterator LhsIterator; + static void run(const Lhs& lhs, Rhs& other) + { + LhsEval lhsEval(lhs); + for(Index col=0 ; col=0; --i) + { + Scalar& tmp = other.coeffRef(i,col); + if (tmp!=Scalar(0)) // optimization when other is actually sparse + { + if(!(Mode & UnitDiag)) + { + // TODO replace this by a binary search. make sure the binary search is safe for partially sorted elements + LhsIterator it(lhsEval, i); + while(it && it.index()!=i) + ++it; + eigen_assert(it && it.index()==i); + other.coeffRef(i,col) /= it.value(); + } + LhsIterator it(lhsEval, i); + for(; it && it.index() +template template -void SparseTriangularView::solveInPlace(MatrixBase& other) const +void TriangularViewImpl::solveInPlace(MatrixBase& other) const { - eigen_assert(m_matrix.cols() == m_matrix.rows() && m_matrix.cols() == other.rows()); + eigen_assert(derived().cols() == derived().rows() && derived().cols() == other.rows()); eigen_assert((!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower))); enum { copy = internal::traits::Flags & RowMajorBit }; @@ -178,21 +332,23 @@ void SparseTriangularView::solveInPlace(MatrixBase::type, OtherDerived&>::type OtherCopy; OtherCopy otherCopy(other.derived()); - internal::sparse_solve_triangular_selector::type, Mode>::run(m_matrix, otherCopy); + internal::sparse_solve_triangular_selector::type, Mode>::run(derived().nestedExpression(), otherCopy); if (copy) other = otherCopy; } -template +#ifndef EIGEN_TEST_EVALUATORS +template template typename internal::plain_matrix_type_column_major::type -SparseTriangularView::solve(const MatrixBase& other) const +TriangularViewImpl::solve(const MatrixBase& other) const { typename internal::plain_matrix_type_column_major::type res(other); solveInPlace(res); return res; } +#endif // pure sparse path @@ -290,11 +446,11 @@ struct sparse_solve_triangular_sparse_selector } // end namespace internal -template +template template -void SparseTriangularView::solveInPlace(SparseMatrixBase& other) const +void TriangularViewImpl::solveInPlace(SparseMatrixBase& other) const { - eigen_assert(m_matrix.cols() == m_matrix.rows() && m_matrix.cols() == other.rows()); + eigen_assert(derived().cols() == derived().rows() && derived().cols() == other.rows()); eigen_assert( (!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower))); // enum { copy = internal::traits::Flags & RowMajorBit }; @@ -303,7 +459,7 @@ void SparseTriangularView::solveInPlace(SparseMatrixBase::type, OtherDerived&>::type OtherCopy; // OtherCopy otherCopy(other.derived()); - internal::sparse_solve_triangular_sparse_selector::run(m_matrix, other.derived()); + internal::sparse_solve_triangular_sparse_selector::run(derived().nestedExpression(), other.derived()); // if (copy) // other = otherCopy;