diff options
31 files changed, 1919 insertions, 199 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index dd74b7674..f304401cb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -76,6 +76,12 @@ if(MSVC) endif(EIGEN_TEST_SSE2) endif(MSVC) +option(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION "Disable explicit vectorization in tests/examples" OFF) +if(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION) + add_definitions(EIGEN_DONT_VECTORIZE=1) + message("Disabling vectorization in tests/examples") +endif(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION) + include_directories(${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR}) add_subdirectory(Eigen) diff --git a/Eigen/Sparse b/Eigen/Sparse index 168511c66..b73946c74 100644 --- a/Eigen/Sparse +++ b/Eigen/Sparse @@ -79,7 +79,15 @@ namespace Eigen { #include "src/Sparse/SparseMatrix.h" #include "src/Sparse/SparseVector.h" #include "src/Sparse/CoreIterators.h" +#include "src/Sparse/SparseTranspose.h" +#include "src/Sparse/SparseCwise.h" +#include "src/Sparse/SparseCwiseUnaryOp.h" +#include "src/Sparse/SparseCwiseBinaryOp.h" +#include "src/Sparse/SparseDot.h" +#include "src/Sparse/SparseAssign.h" #include "src/Sparse/SparseRedux.h" +#include "src/Sparse/SparseFuzzy.h" +#include "src/Sparse/SparseFlagged.h" #include "src/Sparse/SparseProduct.h" #include "src/Sparse/TriangularSolver.h" #include "src/Sparse/SparseLLT.h" diff --git a/Eigen/src/Core/CwiseBinaryOp.h b/Eigen/src/Core/CwiseBinaryOp.h index c5bead928..c4223e220 100644 --- a/Eigen/src/Core/CwiseBinaryOp.h +++ b/Eigen/src/Core/CwiseBinaryOp.h @@ -86,14 +86,12 @@ class CwiseBinaryOp : ei_no_assignment_operator, typedef typename ei_traits<CwiseBinaryOp>::LhsNested LhsNested; typedef typename ei_traits<CwiseBinaryOp>::RhsNested RhsNested; - class InnerIterator; - EIGEN_STRONG_INLINE CwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp()) : m_lhs(lhs), m_rhs(rhs), m_functor(func) { // we require Lhs and Rhs to have the same scalar type. Currently there is no example of a binary functor // that would take two operands of different types. If there were such an example, then this check should be - // moved to the BinaryOp functors, on a per-case basis. This would however require a change in the BinaryOp functors, as + // moved to the BinaryOp functors, on a per-case basis. This would however require a change in the BinaryOp functors, as // currently they take only one typename Scalar template parameter. // It is tempting to always allow mixing different types but remember that this is often impossible in the vectorized paths. // So allowing mixing different types gives very unexpected errors when enabling vectorization, when the user tries to diff --git a/Eigen/src/Core/CwiseUnaryOp.h b/Eigen/src/Core/CwiseUnaryOp.h index 68be2ac54..076d568e0 100644 --- a/Eigen/src/Core/CwiseUnaryOp.h +++ b/Eigen/src/Core/CwiseUnaryOp.h @@ -64,8 +64,6 @@ class CwiseUnaryOp : ei_no_assignment_operator, EIGEN_GENERIC_PUBLIC_INTERFACE(CwiseUnaryOp) - class InnerIterator; - inline CwiseUnaryOp(const MatrixType& mat, const UnaryOp& func = UnaryOp()) : m_matrix(mat), m_functor(func) {} diff --git a/Eigen/src/Core/Dot.h b/Eigen/src/Core/Dot.h index 71e7030dc..b66fcbaae 100644 --- a/Eigen/src/Core/Dot.h +++ b/Eigen/src/Core/Dot.h @@ -143,13 +143,12 @@ struct ei_dot_vec_unroller<Derived1, Derived2, Index, Stop, true> template<typename Derived1, typename Derived2, int Vectorization = ei_dot_traits<Derived1, Derived2>::Vectorization, - int Unrolling = ei_dot_traits<Derived1, Derived2>::Unrolling, - int Storage = (ei_traits<Derived1>::Flags | ei_traits<Derived2>::Flags) & SparseBit + int Unrolling = ei_dot_traits<Derived1, Derived2>::Unrolling > struct ei_dot_impl; template<typename Derived1, typename Derived2> -struct ei_dot_impl<Derived1, Derived2, NoVectorization, NoUnrolling, IsDense> +struct ei_dot_impl<Derived1, Derived2, NoVectorization, NoUnrolling> { typedef typename Derived1::Scalar Scalar; static Scalar run(const Derived1& v1, const Derived2& v2) @@ -164,12 +163,12 @@ struct ei_dot_impl<Derived1, Derived2, NoVectorization, NoUnrolling, IsDense> }; template<typename Derived1, typename Derived2> -struct ei_dot_impl<Derived1, Derived2, NoVectorization, CompleteUnrolling, IsDense> +struct ei_dot_impl<Derived1, Derived2, NoVectorization, CompleteUnrolling> : public ei_dot_novec_unroller<Derived1, Derived2, 0, Derived1::SizeAtCompileTime> {}; template<typename Derived1, typename Derived2> -struct ei_dot_impl<Derived1, Derived2, LinearVectorization, NoUnrolling, IsDense> +struct ei_dot_impl<Derived1, Derived2, LinearVectorization, NoUnrolling> { typedef typename Derived1::Scalar Scalar; typedef typename ei_packet_traits<Scalar>::type PacketScalar; @@ -222,7 +221,7 @@ struct ei_dot_impl<Derived1, Derived2, LinearVectorization, NoUnrolling, IsDense }; template<typename Derived1, typename Derived2> -struct ei_dot_impl<Derived1, Derived2, LinearVectorization, CompleteUnrolling, IsDense> +struct ei_dot_impl<Derived1, Derived2, LinearVectorization, CompleteUnrolling> { typedef typename Derived1::Scalar Scalar; typedef typename ei_packet_traits<Scalar>::type PacketScalar; diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 42d3bece3..eecd24c85 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -251,8 +251,8 @@ template<typename Derived> class MatrixBase { return lazyAssign(other._expression()); } /** Overloaded for sparse product evaluation */ - template<typename Derived1, typename Derived2> - Derived& lazyAssign(const Product<Derived1,Derived2,SparseProduct>& product); + /*template<typename Derived1, typename Derived2> + Derived& lazyAssign(const Product<Derived1,Derived2,SparseProduct>& product);*/ CommaInitializer<Derived> operator<< (const Scalar& s); diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index 2fd9d5875..305efc3dc 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -79,7 +79,6 @@ struct ProductReturnType<Lhs,Rhs,CacheFriendlyProduct> * - NormalProduct * - CacheFriendlyProduct * - DiagonalProduct - * - SparseProduct */ template<typename Lhs, typename Rhs> struct ei_product_mode { @@ -87,8 +86,6 @@ template<typename Lhs, typename Rhs> struct ei_product_mode value = ((Rhs::Flags&Diagonal)==Diagonal) || ((Lhs::Flags&Diagonal)==Diagonal) ? DiagonalProduct - : (Rhs::Flags & Lhs::Flags & SparseBit) - ? SparseProduct : Lhs::MaxColsAtCompileTime == Dynamic && ( Lhs::MaxRowsAtCompileTime == Dynamic || Rhs::MaxColsAtCompileTime == Dynamic ) diff --git a/Eigen/src/Core/Sum.h b/Eigen/src/Core/Sum.h index ab9cdad1f..e30392534 100644 --- a/Eigen/src/Core/Sum.h +++ b/Eigen/src/Core/Sum.h @@ -154,13 +154,12 @@ struct ei_sum_vec_unroller<Derived, Index, Stop, true> template<typename Derived, int Vectorization = ei_sum_traits<Derived>::Vectorization, - int Unrolling = ei_sum_traits<Derived>::Unrolling, - int Storage = ei_traits<Derived>::Flags & SparseBit + int Unrolling = ei_sum_traits<Derived>::Unrolling > struct ei_sum_impl; template<typename Derived> -struct ei_sum_impl<Derived, NoVectorization, NoUnrolling, IsDense> +struct ei_sum_impl<Derived, NoVectorization, NoUnrolling> { typedef typename Derived::Scalar Scalar; static Scalar run(const Derived& mat) @@ -178,12 +177,12 @@ struct ei_sum_impl<Derived, NoVectorization, NoUnrolling, IsDense> }; template<typename Derived> -struct ei_sum_impl<Derived, NoVectorization, CompleteUnrolling, IsDense> +struct ei_sum_impl<Derived, NoVectorization, CompleteUnrolling> : public ei_sum_novec_unroller<Derived, 0, Derived::SizeAtCompileTime> {}; template<typename Derived> -struct ei_sum_impl<Derived, LinearVectorization, NoUnrolling,IsDense> +struct ei_sum_impl<Derived, LinearVectorization, NoUnrolling> { typedef typename Derived::Scalar Scalar; typedef typename ei_packet_traits<Scalar>::type PacketScalar; @@ -228,7 +227,7 @@ struct ei_sum_impl<Derived, LinearVectorization, NoUnrolling,IsDense> }; template<typename Derived> -struct ei_sum_impl<Derived, LinearVectorization, CompleteUnrolling, IsDense> +struct ei_sum_impl<Derived, LinearVectorization, CompleteUnrolling> { typedef typename Derived::Scalar Scalar; typedef typename ei_packet_traits<Scalar>::type PacketScalar; diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h index 1137e42dd..93028ed49 100644 --- a/Eigen/src/Core/Transpose.h +++ b/Eigen/src/Core/Transpose.h @@ -63,8 +63,6 @@ template<typename MatrixType> class Transpose EIGEN_GENERIC_PUBLIC_INTERFACE(Transpose) - class InnerIterator; - inline Transpose(const MatrixType& matrix) : m_matrix(matrix) {} EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Transpose) @@ -185,7 +183,7 @@ struct ei_inplace_transpose_selector<MatrixType,false> { // non square matrix * * In most cases it is probably better to simply use the transposed expression * of a matrix. However, when transposing the matrix data itself is really needed, - * then this "in-place" version is probably the right choice because it provides + * then this "in-place" version is probably the right choice because it provides * the following additional features: * - less error prone: doing the same operation with .transpose() requires special care: * \code m.set(m.transpose().eval()); \endcode diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 10e50f43b..f2c76cc01 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -201,7 +201,7 @@ enum { ForceAligned, AsRequested }; enum { ConditionalJumpCost = 5 }; enum CornerType { TopLeft, TopRight, BottomLeft, BottomRight }; enum DirectionType { Vertical, Horizontal }; -enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, DiagonalProduct, SparseProduct }; +enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, DiagonalProduct }; enum { /** \internal Equivalent to a slice vectorization for fixed-size matrices having good alignment diff --git a/Eigen/src/Sparse/CoreIterators.h b/Eigen/src/Sparse/CoreIterators.h index b1d113a95..d57eac287 100644 --- a/Eigen/src/Sparse/CoreIterators.h +++ b/Eigen/src/Sparse/CoreIterators.h @@ -28,19 +28,26 @@ /* This file contains the respective InnerIterator definition of the expressions defined in Eigen/Core */ -template<typename Derived> -class MatrixBase<Derived>::InnerIterator +/** \class InnerIterator + * \brief An InnerIterator allows to loop over the element of a sparse (or dense) matrix or expression + * + * todo + */ +// template<typename XprType> InnerIterator; + +// generic version for dense matrix and expressions +template<typename Derived> class MatrixBase<Derived>::InnerIterator { typedef typename Derived::Scalar Scalar; public: - EIGEN_STRONG_INLINE InnerIterator(const Derived& mat, int outer) - : m_matrix(mat), m_inner(0), m_outer(outer), m_end(mat.rows()) + EIGEN_STRONG_INLINE InnerIterator(const Derived& expr, int outer) + : m_expression(expr), m_inner(0), m_outer(outer), m_end(expr.rows()) {} EIGEN_STRONG_INLINE Scalar value() const { - return (Derived::Flags&RowMajorBit) ? m_matrix.coeff(m_outer, m_inner) - : m_matrix.coeff(m_inner, m_outer); + return (Derived::Flags&RowMajorBit) ? m_expression.coeff(m_outer, m_inner) + : m_expression.coeff(m_inner, m_outer); } EIGEN_STRONG_INLINE InnerIterator& operator++() { m_inner++; return *this; } @@ -50,22 +57,14 @@ class MatrixBase<Derived>::InnerIterator EIGEN_STRONG_INLINE operator bool() const { return m_inner < m_end && m_inner>=0; } protected: - const Derived& m_matrix; + const Derived& m_expression; int m_inner; const int m_outer; const int m_end; }; -template<typename MatrixType> -class Transpose<MatrixType>::InnerIterator : public MatrixType::InnerIterator -{ - public: - - EIGEN_STRONG_INLINE InnerIterator(const Transpose& trans, int outer) - : MatrixType::InnerIterator(trans.m_matrix, outer) - {} -}; +/* template<typename MatrixType, int BlockRows, int BlockCols, int PacketAccess, int _DirectAccessStatus> class Block<MatrixType, BlockRows, BlockCols, PacketAccess, _DirectAccessStatus>::InnerIterator { @@ -138,50 +137,9 @@ class Block<MatrixType, BlockRows, BlockCols, PacketAccess, IsSparse>::InnerIter int m_start; int m_end; int m_offset; -}; - -template<typename UnaryOp, typename MatrixType> -class CwiseUnaryOp<UnaryOp,MatrixType>::InnerIterator -{ - typedef typename CwiseUnaryOp::Scalar Scalar; - typedef typename ei_traits<CwiseUnaryOp>::_MatrixTypeNested _MatrixTypeNested; - typedef typename _MatrixTypeNested::InnerIterator MatrixTypeIterator; - public: - - EIGEN_STRONG_INLINE InnerIterator(const CwiseUnaryOp& unaryOp, int outer) - : m_iter(unaryOp.m_matrix,outer), m_functor(unaryOp.m_functor), m_id(-1) - { - this->operator++(); - } - - EIGEN_STRONG_INLINE InnerIterator& operator++() - { - if (m_iter) - { - m_id = m_iter.index(); - m_value = m_functor(m_iter.value()); - ++m_iter; - } - else - { - m_id = -1; - } - return *this; - } - - EIGEN_STRONG_INLINE Scalar value() const { return m_value; } - - EIGEN_STRONG_INLINE int index() const { return m_id; } - - EIGEN_STRONG_INLINE operator bool() const { return m_id>=0; } - - protected: - MatrixTypeIterator m_iter; - const UnaryOp& m_functor; - Scalar m_value; - int m_id; -}; +};*/ +/* template<typename T> struct ei_is_scalar_product { enum { ret = false }; }; template<typename T> struct ei_is_scalar_product<ei_scalar_product_op<T> > { enum { ret = true }; }; @@ -204,8 +162,8 @@ class CwiseBinaryOp<BinaryOp,Lhs,Rhs>::InnerIterator EIGEN_STRONG_INLINE InnerIterator(const CwiseBinaryOp& binOp, int outer) : Base(binOp.m_lhs,binOp.m_rhs,binOp.m_functor,outer) {} -}; - +};*/ +/* template<typename BinaryOp, typename Lhs, typename Rhs, typename Derived> class CwiseBinaryOpInnerIterator { @@ -270,7 +228,7 @@ class CwiseBinaryOpInnerIterator Scalar m_value; int m_id; }; -/* + template<typename T, typename Lhs, typename Rhs, typename Derived> class CwiseBinaryOpInnerIterator<ei_scalar_product_op<T>,Lhs,Rhs,Derived> { diff --git a/Eigen/src/Sparse/SparseAssign.h b/Eigen/src/Sparse/SparseAssign.h new file mode 100644 index 000000000..e69de29bb --- /dev/null +++ b/Eigen/src/Sparse/SparseAssign.h diff --git a/Eigen/src/Sparse/SparseBlock.h b/Eigen/src/Sparse/SparseBlock.h index 55590d9e2..48be3754f 100644 --- a/Eigen/src/Sparse/SparseBlock.h +++ b/Eigen/src/Sparse/SparseBlock.h @@ -23,9 +23,76 @@ // License and a copy of the GNU General Public License along with // Eigen. If not, see <http://www.gnu.org/licenses/>. -#ifndef EIGEN_SPARSEBLOCK_H -#define EIGEN_SPARSEBLOCK_H +#ifndef EIGEN_SPARSE_BLOCK_H +#define EIGEN_SPARSE_BLOCK_H +template<typename MatrixType> +struct ei_traits<SparseInnerVector<MatrixType> > +{ + typedef typename ei_traits<MatrixType>::Scalar Scalar; + enum { + IsRowMajor = (int(MatrixType::Flags)&RowMajorBit)==RowMajorBit, + Flags = MatrixType::Flags, + RowsAtCompileTime = IsRowMajor ? 1 : MatrixType::RowsAtCompileTime, + ColsAtCompileTime = IsRowMajor ? MatrixType::ColsAtCompileTime : 1, + CoeffReadCost = MatrixType::CoeffReadCost + }; +}; + +template<typename MatrixType> +class SparseInnerVector : ei_no_assignment_operator, + public SparseMatrixBase<SparseInnerVector<MatrixType> > +{ + enum { + IsRowMajor = ei_traits<SparseInnerVector>::IsRowMajor + }; +public: + + EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseInnerVector) + class InnerIterator; + + inline SparseInnerVector(const MatrixType& matrix, int outer) + : m_matrix(matrix), m_outer(outer) + { + ei_assert( (outer>=0) && (outer<matrix.outerSize()) ); + } + + EIGEN_STRONG_INLINE int rows() const { return IsRowMajor ? 1 : m_matrix.rows(); } + EIGEN_STRONG_INLINE int cols() const { return IsRowMajor ? m_matrix.cols() : 1; } + + protected: + + const typename MatrixType::Nested m_matrix; + int m_outer; + +}; + +template<typename MatrixType> +class SparseInnerVector<MatrixType>::InnerIterator : public MatrixType::InnerIterator +{ +public: + inline InnerIterator(const SparseInnerVector& xpr, int outer=0) + : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outer) + { + ei_assert(outer==0); + } +}; + +/** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this + * is col-major (resp. row-major). + */ +template<typename Derived> +SparseInnerVector<Derived> SparseMatrixBase<Derived>::innerVector(int outer) +{ return SparseInnerVector<Derived>(derived(), outer); } + +/** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this + * is col-major (resp. row-major). Read-only. + */ +template<typename Derived> +const SparseInnerVector<Derived> SparseMatrixBase<Derived>::innerVector(int outer) const +{ return SparseInnerVector<Derived>(derived(), outer); } + +# if 0 template<typename MatrixType, int BlockRows, int BlockCols, int PacketAccess> class Block<MatrixType,BlockRows,BlockCols,PacketAccess,IsSparse> : public SparseMatrixBase<Block<MatrixType,BlockRows,BlockCols,PacketAccess,IsSparse> > @@ -117,6 +184,6 @@ public: const ei_int_if_dynamic<ColsAtCompileTime> m_blockCols; }; +#endif - -#endif // EIGEN_SPARSEBLOCK_H +#endif // EIGEN_SPARSE_BLOCK_H diff --git a/Eigen/src/Sparse/SparseCwise.h b/Eigen/src/Sparse/SparseCwise.h new file mode 100644 index 000000000..0697439d3 --- /dev/null +++ b/Eigen/src/Sparse/SparseCwise.h @@ -0,0 +1,181 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_SPARSE_CWISE_H +#define EIGEN_SPARSE_CWISE_H + +/** \internal + * convenient macro to defined the return type of a cwise binary operation */ +#define EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(OP) \ + CwiseBinaryOp<OP<typename ei_traits<ExpressionType>::Scalar>, ExpressionType, OtherDerived> + +#define EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE \ + SparseCwiseBinaryOp< \ + ei_scalar_product_op< \ + typename ei_scalar_product_traits< \ + typename ei_traits<ExpressionType>::Scalar, \ + typename ei_traits<OtherDerived>::Scalar \ + >::ReturnType \ + >, \ + ExpressionType, \ + OtherDerived \ + > + +/** \internal + * convenient macro to defined the return type of a cwise unary operation */ +#define EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(OP) \ + SparseCwiseUnaryOp<OP<typename ei_traits<ExpressionType>::Scalar>, ExpressionType> + +/** \internal + * convenient macro to defined the return type of a cwise comparison to a scalar */ +/*#define EIGEN_CWISE_COMP_TO_SCALAR_RETURN_TYPE(OP) \ + CwiseBinaryOp<OP<typename ei_traits<ExpressionType>::Scalar>, ExpressionType, \ + NestByValue<typename ExpressionType::ConstantReturnType> >*/ + +template<typename ExpressionType> class SparseCwise +{ + public: + + typedef typename ei_traits<ExpressionType>::Scalar Scalar; + typedef typename ei_meta_if<ei_must_nest_by_value<ExpressionType>::ret, + ExpressionType, const ExpressionType&>::ret ExpressionTypeNested; + typedef CwiseUnaryOp<ei_scalar_add_op<Scalar>, ExpressionType> ScalarAddReturnType; + + inline SparseCwise(const ExpressionType& matrix) : m_matrix(matrix) {} + + /** \internal */ + inline const ExpressionType& _expression() const { return m_matrix; } + + template<typename OtherDerived> + const EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE + operator*(const SparseMatrixBase<OtherDerived> &other) const; + + template<typename OtherDerived> + const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op) + operator/(const SparseMatrixBase<OtherDerived> &other) const; + + template<typename OtherDerived> + const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_min_op) + min(const SparseMatrixBase<OtherDerived> &other) const; + + template<typename OtherDerived> + const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_max_op) + max(const SparseMatrixBase<OtherDerived> &other) const; + + const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_abs_op) abs() const; + const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_abs2_op) abs2() const; +// const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_square_op) square() const; +// const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_cube_op) cube() const; +// const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_inverse_op) inverse() const; +// const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_sqrt_op) sqrt() const; +// const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_exp_op) exp() const; +// const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_log_op) log() const; +// const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_cos_op) cos() const; +// const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_sin_op) sin() const; +// const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_pow_op) pow(const Scalar& exponent) const; + + /* + const ScalarAddReturnType + operator+(const Scalar& scalar) const; + + friend const ScalarAddReturnType + operator+(const Scalar& scalar, const Cwise& mat) + { return mat + scalar; } + + ExpressionType& operator+=(const Scalar& scalar); + + const ScalarAddReturnType + operator-(const Scalar& scalar) const; + + ExpressionType& operator-=(const Scalar& scalar); + + template<typename OtherDerived> + inline ExpressionType& operator*=(const MatrixBase<OtherDerived> &other); + + template<typename OtherDerived> + inline ExpressionType& operator/=(const MatrixBase<OtherDerived> &other); + + template<typename OtherDerived> const EIGEN_CWISE_BINOP_RETURN_TYPE(std::less) + operator<(const MatrixBase<OtherDerived>& other) const; + + template<typename OtherDerived> const EIGEN_CWISE_BINOP_RETURN_TYPE(std::less_equal) + operator<=(const MatrixBase<OtherDerived>& other) const; + + template<typename OtherDerived> const EIGEN_CWISE_BINOP_RETURN_TYPE(std::greater) + operator>(const MatrixBase<OtherDerived>& other) const; + + template<typename OtherDerived> const EIGEN_CWISE_BINOP_RETURN_TYPE(std::greater_equal) + operator>=(const MatrixBase<OtherDerived>& other) const; + + template<typename OtherDerived> const EIGEN_CWISE_BINOP_RETURN_TYPE(std::equal_to) + operator==(const MatrixBase<OtherDerived>& other) const; + + template<typename OtherDerived> const EIGEN_CWISE_BINOP_RETURN_TYPE(std::not_equal_to) + operator!=(const MatrixBase<OtherDerived>& other) const; + + // comparisons to a scalar value + const EIGEN_CWISE_COMP_TO_SCALAR_RETURN_TYPE(std::less) + operator<(Scalar s) const; + + const EIGEN_CWISE_COMP_TO_SCALAR_RETURN_TYPE(std::less_equal) + operator<=(Scalar s) const; + + const EIGEN_CWISE_COMP_TO_SCALAR_RETURN_TYPE(std::greater) + operator>(Scalar s) const; + + const EIGEN_CWISE_COMP_TO_SCALAR_RETURN_TYPE(std::greater_equal) + operator>=(Scalar s) const; + + const EIGEN_CWISE_COMP_TO_SCALAR_RETURN_TYPE(std::equal_to) + operator==(Scalar s) const; + + const EIGEN_CWISE_COMP_TO_SCALAR_RETURN_TYPE(std::not_equal_to) + operator!=(Scalar s) const; + */ + + // allow to extend SparseCwise outside Eigen + #ifdef EIGEN_SPARSE_CWISE_PLUGIN + #include EIGEN_SPARSE_CWISE_PLUGIN + #endif + + protected: + ExpressionTypeNested m_matrix; +}; + +template<typename Derived> +inline const SparseCwise<Derived> +SparseMatrixBase<Derived>::cwise() const +{ + return derived(); +} + +template<typename Derived> +inline SparseCwise<Derived> +SparseMatrixBase<Derived>::cwise() +{ + return derived(); +} + +#endif // EIGEN_SPARSE_CWISE_H diff --git a/Eigen/src/Sparse/SparseCwiseBinaryOp.h b/Eigen/src/Sparse/SparseCwiseBinaryOp.h new file mode 100644 index 000000000..99bbae0cd --- /dev/null +++ b/Eigen/src/Sparse/SparseCwiseBinaryOp.h @@ -0,0 +1,560 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_SPARSE_CWISE_BINARY_OP_H +#define EIGEN_SPARSE_CWISE_BINARY_OP_H + +// 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 + +// template<typename BinaryOp, typename Lhs, typename Rhs> class ei_sparse_cwise_binary_op +// : ei_no_assignment_operator +// { +// typedef CwiseBinaryOp<BinaryOp,Lhs,Rhs> CwiseBinaryXpr; +// public: +// +// typedef typename ei_traits<CwiseBinaryXpr>::LhsNested LhsNested; +// typedef typename ei_traits<CwiseBinaryXpr>::RhsNested RhsNested; +// +// EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp()) +// : m_lhs(lhs), m_rhs(rhs), m_functor(func) +// { +// EIGEN_STATIC_ASSERT((ei_functor_allows_mixing_real_and_complex<BinaryOp>::ret +// ? int(ei_is_same_type<typename Lhs::RealScalar, typename Rhs::RealScalar>::ret) +// : int(ei_is_same_type<typename Lhs::Scalar, typename Rhs::Scalar>::ret)), +// YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) +// // require the sizes to match +// EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Lhs, Rhs) +// ei_assert(lhs.rows() == rhs.rows() && lhs.cols() == rhs.cols()); +// } +// +// EIGEN_STRONG_INLINE int rows() const { return m_lhs.rows(); } +// EIGEN_STRONG_INLINE int cols() const { return m_lhs.cols(); } +// +// EIGEN_STRONG_INLINE const LhsNested& _lhs() const { return m_lhs; } +// EIGEN_STRONG_INLINE const RhsNested& _rhs() const { return m_rhs; } +// EIGEN_STRONG_INLINE const UnaryOp& _functor() const { return m_functor; } +// +// protected: +// const LhsNested m_lhs; +// const RhsNested m_rhs; +// const BinaryOp m_functor; +// } + +template<typename BinaryOp, typename Lhs, typename Rhs> +struct ei_traits<SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs> > +{ + typedef typename ei_result_of< + BinaryOp( + typename Lhs::Scalar, + typename Rhs::Scalar + ) + >::type Scalar; + typedef typename Lhs::Nested LhsNested; + typedef typename Rhs::Nested RhsNested; + typedef typename ei_unref<LhsNested>::type _LhsNested; + typedef typename ei_unref<RhsNested>::type _RhsNested; + enum { + LhsCoeffReadCost = _LhsNested::CoeffReadCost, + RhsCoeffReadCost = _RhsNested::CoeffReadCost, + LhsFlags = _LhsNested::Flags, + RhsFlags = _RhsNested::Flags, + RowsAtCompileTime = Lhs::RowsAtCompileTime, + ColsAtCompileTime = Lhs::ColsAtCompileTime, + MaxRowsAtCompileTime = Lhs::MaxRowsAtCompileTime, + MaxColsAtCompileTime = Lhs::MaxColsAtCompileTime, + Flags = (int(LhsFlags) | int(RhsFlags)) & HereditaryBits, + CoeffReadCost = LhsCoeffReadCost + RhsCoeffReadCost + ei_functor_traits<BinaryOp>::Cost + }; +}; + +template<typename BinaryOp, typename Lhs, typename Rhs> +class SparseCwiseBinaryOp : ei_no_assignment_operator, + public SparseMatrixBase<SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs> > +{ + public: + + class InnerIterator; + + EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseCwiseBinaryOp) + typedef typename ei_traits<SparseCwiseBinaryOp>::LhsNested LhsNested; + typedef typename ei_traits<SparseCwiseBinaryOp>::RhsNested RhsNested; + typedef typename ei_unref<LhsNested>::type _LhsNested; + typedef typename ei_unref<RhsNested>::type _RhsNested; + + EIGEN_STRONG_INLINE SparseCwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp()) + : m_lhs(lhs), m_rhs(rhs), m_functor(func) + { + EIGEN_STATIC_ASSERT((ei_functor_allows_mixing_real_and_complex<BinaryOp>::ret + ? int(ei_is_same_type<typename Lhs::RealScalar, typename Rhs::RealScalar>::ret) + : int(ei_is_same_type<typename Lhs::Scalar, typename Rhs::Scalar>::ret)), + YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) + // require the sizes to match + EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Lhs, Rhs) + ei_assert(lhs.rows() == rhs.rows() && lhs.cols() == rhs.cols()); + } + + EIGEN_STRONG_INLINE int rows() const { return m_lhs.rows(); } + EIGEN_STRONG_INLINE int cols() const { return m_lhs.cols(); } + + EIGEN_STRONG_INLINE const _LhsNested& _lhs() const { return m_lhs; } + EIGEN_STRONG_INLINE const _RhsNested& _rhs() const { return m_rhs; } + EIGEN_STRONG_INLINE const BinaryOp& _functor() const { return m_functor; } + + protected: + const LhsNested m_lhs; + const RhsNested m_rhs; + const BinaryOp m_functor; +}; +/* +template<typename BinaryOp, typename Lhs, typename Rhs> +class CwiseBinaryOp<BinaryOp, SparseMatrixBase<Lhs>, Rhs> + : public ei_sparse_cwise_binary_op<BinaryOp, Lhs, Rhs>, + public SparseMatrixBase<CwiseBinaryOp<BinaryOp, SparseMatrixBase<Lhs>, Rhs> > +{ + public: + EIGEN_GENERIC_PUBLIC_INTERFACE(CwiseBinaryOp) + EIGEN_STRONG_INLINE CwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp()) + : ei_sparse_cwise_binary_op<BinaryOp, Lhs, Rhs>(lhs, rhs, func) + {} +}; +template<typename BinaryOp, typename Lhs, typename Rhs> +class CwiseBinaryOp<BinaryOp, SparseMatrixBase<Lhs>, SparseMatrixBase<Rhs> > + : public ei_sparse_cwise_binary_op<BinaryOp, Lhs, Rhs>, + public SparseMatrixBase<CwiseBinaryOp<BinaryOp, Lhs, Rhs> > +{ + public: + EIGEN_GENERIC_PUBLIC_INTERFACE(CwiseBinaryOp) + EIGEN_STRONG_INLINE CwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp()) + : ei_sparse_cwise_binary_op<BinaryOp, Lhs, Rhs>(lhs, rhs, func) + {} +};*/ + +// template<typename BinaryOp, typename Lhs, typename RhsDerived> +// class CwiseBinaryOp<BinaryOp, Lhs, SparseMatrixBase<Rhs> > : ei_no_assignment_operator, +// public SparseMatrixBase<CwiseBinaryOp<BinaryOp, SparseMatrixBase<LhsDerived>, SparseMatrixBase<Rhs> > > +// { +// public: +// +// EIGEN_GENERIC_PUBLIC_INTERFACE(CwiseBinaryOp) +// typedef typename ei_traits<CwiseBinaryOp>::LhsNested LhsNested; +// typedef typename ei_traits<CwiseBinaryOp>::RhsNested RhsNested; +// +// EIGEN_STRONG_INLINE CwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp()) +// : m_lhs(lhs), m_rhs(rhs), m_functor(func) +// { +// EIGEN_STATIC_ASSERT((ei_functor_allows_mixing_real_and_complex<BinaryOp>::ret +// ? int(ei_is_same_type<typename Lhs::RealScalar, typename Rhs::RealScalar>::ret) +// : int(ei_is_same_type<typename Lhs::Scalar, typename Rhs::Scalar>::ret)), +// YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) +// // require the sizes to match +// EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Lhs, Rhs) +// ei_assert(lhs.rows() == rhs.rows() && lhs.cols() == rhs.cols()); +// } +// +// EIGEN_STRONG_INLINE int rows() const { return m_lhs.rows(); } +// EIGEN_STRONG_INLINE int cols() const { return m_lhs.cols(); } +// +// protected: +// const LhsNested m_lhs; +// const RhsNested m_rhs; +// const BinaryOp m_functor; +// }; +// +// template<typename BinaryOp, typename LhsDerived, typename RhsDerived> +// class CwiseBinaryOp<BinaryOp, SparseMatrixBase<LhsDerived>, SparseMatrixBase<Rhs> > : ei_no_assignment_operator, +// public SparseMatrixBase<CwiseBinaryOp<BinaryOp, SparseMatrixBase<LhsDerived>, SparseMatrixBase<Rhs> > > +// { +// public: +// +// EIGEN_GENERIC_PUBLIC_INTERFACE(CwiseBinaryOp) +// typedef typename ei_traits<CwiseBinaryOp>::LhsNested LhsNested; +// typedef typename ei_traits<CwiseBinaryOp>::RhsNested RhsNested; +// +// EIGEN_STRONG_INLINE CwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp()) +// : m_lhs(lhs), m_rhs(rhs), m_functor(func) +// { +// EIGEN_STATIC_ASSERT((ei_functor_allows_mixing_real_and_complex<BinaryOp>::ret +// ? int(ei_is_same_type<typename Lhs::RealScalar, typename Rhs::RealScalar>::ret) +// : int(ei_is_same_type<typename Lhs::Scalar, typename Rhs::Scalar>::ret)), +// YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) +// // require the sizes to match +// EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Lhs, Rhs) +// ei_assert(lhs.rows() == rhs.rows() && lhs.cols() == rhs.cols()); +// } +// +// EIGEN_STRONG_INLINE int rows() const { return m_lhs.rows(); } +// EIGEN_STRONG_INLINE int cols() const { return m_lhs.cols(); } +// +// protected: +// const LhsNested m_lhs; +// const RhsNested m_rhs; +// const BinaryOp m_functor; +// }; + +// template<typename T> struct ei_is_scalar_product { enum { ret = false }; }; +// template<typename T> struct ei_is_scalar_product<ei_scalar_product_op<T> > { enum { ret = true }; }; + +template<typename BinaryOp, typename Lhs, typename Rhs, typename Derived, + int _LhsStorageMode = int(Lhs::Flags) & SparseBit, + int _RhsStorageMode = int(Rhs::Flags) & SparseBit> +class ei_sparse_cwise_binary_op_inner_iterator_selector; + +template<typename BinaryOp, typename Lhs, typename Rhs> +class SparseCwiseBinaryOp<BinaryOp,Lhs,Rhs>::InnerIterator + : public ei_sparse_cwise_binary_op_inner_iterator_selector<BinaryOp,Lhs,Rhs, typename SparseCwiseBinaryOp<BinaryOp,Lhs,Rhs>::InnerIterator> +{ + typedef ei_sparse_cwise_binary_op_inner_iterator_selector< + BinaryOp,Lhs,Rhs, typename SparseCwiseBinaryOp<BinaryOp,Lhs,Rhs>::InnerIterator> Base; + public: + typedef typename SparseCwiseBinaryOp::Scalar Scalar; + typedef typename ei_traits<SparseCwiseBinaryOp>::_LhsNested _LhsNested; + typedef typename _LhsNested::InnerIterator LhsIterator; + typedef typename ei_traits<SparseCwiseBinaryOp>::_RhsNested _RhsNested; + typedef typename _RhsNested::InnerIterator RhsIterator; + + EIGEN_STRONG_INLINE InnerIterator(const SparseCwiseBinaryOp& binOp, int outer) + : Base(binOp,outer) + {} +}; + +/*************************************************************************** +* Implementation of inner-iterators +***************************************************************************/ + +// template<typename T> struct ei_is_scalar_product { enum { ret = false }; }; +// template<typename T> struct ei_is_scalar_product<ei_scalar_product_op<T> > { enum { ret = true }; }; + +// helper class + + +// sparse - sparse (generic) +template<typename BinaryOp, typename Lhs, typename Rhs, typename Derived> +class ei_sparse_cwise_binary_op_inner_iterator_selector<BinaryOp, Lhs, Rhs, Derived, IsSparse, IsSparse> +{ + typedef SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs> CwiseBinaryXpr; + typedef typename ei_traits<CwiseBinaryXpr>::Scalar Scalar; + typedef typename ei_traits<CwiseBinaryXpr>::_LhsNested _LhsNested; + typedef typename ei_traits<CwiseBinaryXpr>::_RhsNested _RhsNested; + typedef typename _LhsNested::InnerIterator LhsIterator; + typedef typename _RhsNested::InnerIterator RhsIterator; + public: + + EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, int outer) + : m_lhsIter(xpr._lhs(),outer), m_rhsIter(xpr._rhs()), m_functor(xpr._functor()) + { + this->operator++(); + } + + EIGEN_STRONG_INLINE Derived& 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_id = -1; + } + return *static_cast<Derived*>(this); + } + + EIGEN_STRONG_INLINE Scalar value() const { return m_value; } + + EIGEN_STRONG_INLINE int index() const { return m_id; } + + EIGEN_STRONG_INLINE operator bool() const { return m_id>=0; } + + protected: + LhsIterator m_lhsIter; + RhsIterator m_rhsIter; + const BinaryOp& m_functor; + Scalar m_value; + int m_id; +}; + +// template<typename BinaryOp, typename Lhs, typename Rhs> +// class SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs>::InnerIterator +// : public ei_sparse_cwise_binary_op_inner_iterator_selector<BinaryOp,Lhs,Rhs, typename SparseCwiseBinaryOp<BinaryOp,Lhs,Rhs>::InnerIterator> +// { +// typedef CwiseBinaryOpInnerIterator< +// BinaryOp,Lhs,Rhs, typename CwiseBinaryOp<BinaryOp,Lhs,Rhs>::InnerIterator> Base; +// public: +// typedef typename CwiseBinaryOp::Scalar Scalar; +// typedef typename ei_traits<CwiseBinaryOp>::_LhsNested _LhsNested; +// typedef typename _LhsNested::InnerIterator LhsIterator; +// typedef typename ei_traits<CwiseBinaryOp>::_RhsNested _RhsNested; +// typedef typename _RhsNested::InnerIterator RhsIterator; +// // public: +// EIGEN_STRONG_INLINE InnerIterator(const CwiseBinaryOp& binOp, int outer) +// : Base(binOp.m_lhs,binOp.m_rhs,binOp.m_functor,outer) +// {} +// }; + +// sparse - sparse (product) +template<typename T, typename Lhs, typename Rhs, typename Derived> +class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>, Lhs, Rhs, Derived, IsSparse, IsSparse> +{ + typedef ei_scalar_product_op<T> BinaryOp; + typedef SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs> CwiseBinaryXpr; + typedef typename CwiseBinaryXpr::Scalar Scalar; + typedef typename ei_traits<CwiseBinaryXpr>::_LhsNested _LhsNested; + typedef typename _LhsNested::InnerIterator LhsIterator; + typedef typename ei_traits<CwiseBinaryXpr>::_RhsNested _RhsNested; + typedef typename _RhsNested::InnerIterator RhsIterator; + public: + + EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, int outer) + : m_lhsIter(xpr._lhs(),outer), m_rhsIter(xpr._rhs()), m_functor(xpr._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 Derived& operator++() + { + while (m_lhsIter && m_rhsIter) + { + if (m_lhsIter.index() == m_rhsIter.index()) + { + ++m_lhsIter; + ++m_rhsIter; + break; + } + else if (m_lhsIter.index() < m_rhsIter.index()) + ++m_lhsIter; + else + ++m_rhsIter; + } + return *static_cast<Derived*>(this); + } + + EIGEN_STRONG_INLINE Scalar value() const { return m_functor(m_lhsIter.value(), m_rhsIter.value()); } + + EIGEN_STRONG_INLINE int index() const { return m_lhsIter.index(); } + + EIGEN_STRONG_INLINE operator bool() const { return m_lhsIter && m_rhsIter; } + + protected: + LhsIterator m_lhsIter; + RhsIterator m_rhsIter; + const BinaryOp& m_functor; +}; + +// sparse - dense (product) +template<typename T, typename Lhs, typename Rhs, typename Derived> +class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>, Lhs, Rhs, Derived, IsSparse, IsDense> +{ + typedef ei_scalar_product_op<T> BinaryOp; + typedef SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs> CwiseBinaryXpr; + typedef typename CwiseBinaryXpr::Scalar Scalar; + typedef typename ei_traits<CwiseBinaryXpr>::_LhsNested _LhsNested; + typedef typename _LhsNested::InnerIterator LhsIterator; + enum { IsRowMajor = (int(Lhs::Flags)&RowMajorBit)==RowMajorBit }; + public: + + EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, int outer) + : m_xpr(xpr), m_lhsIter(xpr.lhs(),outer), m_functor(xpr.functor()), m_outer(outer) + {} + + EIGEN_STRONG_INLINE Derived& operator++() + { + ++m_lhsIter; + return *static_cast<Derived*>(this); + } + + EIGEN_STRONG_INLINE Scalar value() const + { return m_functor(m_lhsIter.value(), m_xpr.rhs().coeff(IsRowMajor?m_outer:m_lhsIter.index(),IsRowMajor?m_lhsIter.index():m_outer)); } + + EIGEN_STRONG_INLINE int index() const { return m_lhsIter.index(); } + + EIGEN_STRONG_INLINE operator bool() const { return m_lhsIter; } + + protected: + const CwiseBinaryXpr& m_xpr; + LhsIterator m_lhsIter; + const BinaryOp& m_functor; + const int m_outer; +}; + +// sparse - dense (product) +template<typename T, typename Lhs, typename Rhs, typename Derived> +class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>, Lhs, Rhs, Derived, IsDense, IsSparse> +{ + typedef ei_scalar_product_op<T> BinaryOp; + typedef SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs> CwiseBinaryXpr; + typedef typename CwiseBinaryXpr::Scalar Scalar; + typedef typename ei_traits<CwiseBinaryXpr>::_RhsNested _RhsNested; + typedef typename _RhsNested::InnerIterator RhsIterator; + enum { IsRowMajor = (int(Rhs::Flags)&RowMajorBit)==RowMajorBit }; + public: + + EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, int outer) + : m_xpr(xpr), m_rhsIter(xpr.rhs(),outer), m_functor(xpr.functor()), m_outer(outer) + {} + + EIGEN_STRONG_INLINE Derived& operator++() + { + ++m_rhsIter; + return *static_cast<Derived*>(this); + } + + EIGEN_STRONG_INLINE Scalar value() const + { return m_functor(m_xpr.lhs().coeff(IsRowMajor?m_outer:m_rhsIter.index(),IsRowMajor?m_rhsIter.index():m_outer), m_rhsIter.value()); } + + EIGEN_STRONG_INLINE int index() const { return m_rhsIter.index(); } + + EIGEN_STRONG_INLINE operator bool() const { return m_rhsIter; } + + protected: + const CwiseBinaryXpr& m_xpr; + RhsIterator m_rhsIter; + const BinaryOp& m_functor; + const int m_outer; +}; + + +template<typename Derived> +template<typename OtherDerived> +EIGEN_STRONG_INLINE const SparseCwiseBinaryOp<ei_scalar_difference_op<typename ei_traits<Derived>::Scalar>, + Derived, OtherDerived> +SparseMatrixBase<Derived>::operator-(const SparseMatrixBase<OtherDerived> &other) const +{ + return SparseCwiseBinaryOp<ei_scalar_difference_op<Scalar>, + Derived, OtherDerived>(derived(), other.derived()); +} + +template<typename Derived> +template<typename OtherDerived> +EIGEN_STRONG_INLINE Derived & +SparseMatrixBase<Derived>::operator-=(const SparseMatrixBase<OtherDerived> &other) +{ + return *this = *this - other; +} + +template<typename Derived> +template<typename OtherDerived> +EIGEN_STRONG_INLINE const SparseCwiseBinaryOp<ei_scalar_sum_op<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived> +SparseMatrixBase<Derived>::operator+(const SparseMatrixBase<OtherDerived> &other) const +{ + return SparseCwiseBinaryOp<ei_scalar_sum_op<Scalar>, Derived, OtherDerived>(derived(), other.derived()); +} + +template<typename Derived> +template<typename OtherDerived> +EIGEN_STRONG_INLINE Derived & +SparseMatrixBase<Derived>::operator+=(const SparseMatrixBase<OtherDerived>& other) +{ + return *this = *this + other; +} + +template<typename ExpressionType> +template<typename OtherDerived> +EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE +SparseCwise<ExpressionType>::operator*(const SparseMatrixBase<OtherDerived> &other) const +{ + return EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE(_expression(), other.derived()); +} + +template<typename ExpressionType> +template<typename OtherDerived> +EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op) +SparseCwise<ExpressionType>::operator/(const SparseMatrixBase<OtherDerived> &other) const +{ + return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op)(_expression(), other.derived()); +} + +// template<typename ExpressionType> +// template<typename OtherDerived> +// inline ExpressionType& Cwise<ExpressionType>::operator*=(const SparseMatrixBase<OtherDerived> &other) +// { +// return m_matrix.const_cast_derived() = *this * other; +// } + +// template<typename ExpressionType> +// template<typename OtherDerived> +// inline ExpressionType& Cwise<ExpressionType>::operator/=(const SparseMatrixBase<OtherDerived> &other) +// { +// return m_matrix.const_cast_derived() = *this / other; +// } + +template<typename ExpressionType> +template<typename OtherDerived> +EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_min_op) +SparseCwise<ExpressionType>::min(const SparseMatrixBase<OtherDerived> &other) const +{ + return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_min_op)(_expression(), other.derived()); +} + +template<typename ExpressionType> +template<typename OtherDerived> +EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_max_op) +SparseCwise<ExpressionType>::max(const SparseMatrixBase<OtherDerived> &other) const +{ + return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_max_op)(_expression(), other.derived()); +} + +// template<typename Derived> +// template<typename CustomBinaryOp, typename OtherDerived> +// EIGEN_STRONG_INLINE const CwiseBinaryOp<CustomBinaryOp, Derived, OtherDerived> +// SparseMatrixBase<Derived>::binaryExpr(const SparseMatrixBase<OtherDerived> &other, const CustomBinaryOp& func) const +// { +// return CwiseBinaryOp<CustomBinaryOp, Derived, OtherDerived>(derived(), other.derived(), func); +// } + +#endif // EIGEN_SPARSE_CWISE_BINARY_OP_H diff --git a/Eigen/src/Sparse/SparseCwiseUnaryOp.h b/Eigen/src/Sparse/SparseCwiseUnaryOp.h new file mode 100644 index 000000000..2b2050edd --- /dev/null +++ b/Eigen/src/Sparse/SparseCwiseUnaryOp.h @@ -0,0 +1,175 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_SPARSE_CWISE_UNARY_OP_H +#define EIGEN_SPARSE_CWISE_UNARY_OP_H + +template<typename UnaryOp, typename MatrixType> +struct ei_traits<SparseCwiseUnaryOp<UnaryOp, MatrixType> > : ei_traits<MatrixType> +{ + typedef typename ei_result_of< + UnaryOp(typename MatrixType::Scalar) + >::type Scalar; + typedef typename MatrixType::Nested MatrixTypeNested; + typedef typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested; + enum { + CoeffReadCost = _MatrixTypeNested::CoeffReadCost + ei_functor_traits<UnaryOp>::Cost + }; +}; + +template<typename UnaryOp, typename MatrixType> +class SparseCwiseUnaryOp : ei_no_assignment_operator, + public SparseMatrixBase<SparseCwiseUnaryOp<UnaryOp, MatrixType> > +{ + public: + + class InnerIterator; +// typedef typename ei_unref<LhsNested>::type _LhsNested; + + EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseCwiseUnaryOp) + + inline SparseCwiseUnaryOp(const MatrixType& mat, const UnaryOp& func = UnaryOp()) + : m_matrix(mat), m_functor(func) {} + + EIGEN_STRONG_INLINE int rows() const { return m_matrix.rows(); } + EIGEN_STRONG_INLINE int cols() const { return m_matrix.cols(); } + +// EIGEN_STRONG_INLINE const typename MatrixType::Nested& _matrix() const { return m_matrix; } +// EIGEN_STRONG_INLINE const UnaryOp& _functor() const { return m_functor; } + + protected: + const typename MatrixType::Nested m_matrix; + const UnaryOp m_functor; +}; + + +template<typename UnaryOp, typename MatrixType> +class SparseCwiseUnaryOp<UnaryOp,MatrixType>::InnerIterator +{ + typedef typename SparseCwiseUnaryOp::Scalar Scalar; + typedef typename ei_traits<SparseCwiseUnaryOp>::_MatrixTypeNested _MatrixTypeNested; + typedef typename _MatrixTypeNested::InnerIterator MatrixTypeIterator; + public: + + EIGEN_STRONG_INLINE InnerIterator(const SparseCwiseUnaryOp& unaryOp, int outer) + : m_iter(unaryOp.m_matrix,outer), m_functor(unaryOp.m_functor) + {} + + EIGEN_STRONG_INLINE InnerIterator& operator++() + { ++m_iter; return *this; } + + EIGEN_STRONG_INLINE Scalar value() const { return m_functor(m_iter.value()); } + + EIGEN_STRONG_INLINE int index() const { return m_iter.index(); } + + EIGEN_STRONG_INLINE operator bool() const { return m_iter; } + + protected: + MatrixTypeIterator m_iter; + const UnaryOp& m_functor; +}; + +template<typename Derived> +template<typename CustomUnaryOp> +EIGEN_STRONG_INLINE const SparseCwiseUnaryOp<CustomUnaryOp, Derived> +SparseMatrixBase<Derived>::unaryExpr(const CustomUnaryOp& func) const +{ + return SparseCwiseUnaryOp<CustomUnaryOp, Derived>(derived(), func); +} + +template<typename Derived> +EIGEN_STRONG_INLINE const SparseCwiseUnaryOp<ei_scalar_opposite_op<typename ei_traits<Derived>::Scalar>,Derived> +SparseMatrixBase<Derived>::operator-() const +{ + return derived(); +} + +template<typename ExpressionType> +EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_abs_op) +SparseCwise<ExpressionType>::abs() const +{ + return _expression(); +} + +template<typename ExpressionType> +EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_abs2_op) +SparseCwise<ExpressionType>::abs2() const +{ + return _expression(); +} + +template<typename Derived> +EIGEN_STRONG_INLINE typename SparseMatrixBase<Derived>::ConjugateReturnType +SparseMatrixBase<Derived>::conjugate() const +{ + return ConjugateReturnType(derived()); +} + +template<typename Derived> +EIGEN_STRONG_INLINE const typename SparseMatrixBase<Derived>::RealReturnType +SparseMatrixBase<Derived>::real() const { return derived(); } + +template<typename Derived> +EIGEN_STRONG_INLINE const typename SparseMatrixBase<Derived>::ImagReturnType +SparseMatrixBase<Derived>::imag() const { return derived(); } + +template<typename Derived> +template<typename NewType> +EIGEN_STRONG_INLINE const SparseCwiseUnaryOp<ei_scalar_cast_op<typename ei_traits<Derived>::Scalar, NewType>, Derived> +SparseMatrixBase<Derived>::cast() const +{ + return derived(); +} + +template<typename Derived> +EIGEN_STRONG_INLINE const SparseCwiseUnaryOp<ei_scalar_multiple_op<typename ei_traits<Derived>::Scalar>, Derived> +SparseMatrixBase<Derived>::operator*(const Scalar& scalar) const +{ + return SparseCwiseUnaryOp<ei_scalar_multiple_op<Scalar>, Derived> + (derived(), ei_scalar_multiple_op<Scalar>(scalar)); +} + +template<typename Derived> +EIGEN_STRONG_INLINE const SparseCwiseUnaryOp<ei_scalar_quotient1_op<typename ei_traits<Derived>::Scalar>, Derived> +SparseMatrixBase<Derived>::operator/(const Scalar& scalar) const +{ + return SparseCwiseUnaryOp<ei_scalar_quotient1_op<Scalar>, Derived> + (derived(), ei_scalar_quotient1_op<Scalar>(scalar)); +} + +template<typename Derived> +EIGEN_STRONG_INLINE Derived& +SparseMatrixBase<Derived>::operator*=(const Scalar& other) +{ + return *this = *this * other; +} + +template<typename Derived> +EIGEN_STRONG_INLINE Derived& +SparseMatrixBase<Derived>::operator/=(const Scalar& other) +{ + return *this = *this / other; +} + +#endif // EIGEN_SPARSE_CWISE_UNARY_OP_H diff --git a/Eigen/src/Sparse/SparseDot.h b/Eigen/src/Sparse/SparseDot.h new file mode 100644 index 000000000..7a26e0f4b --- /dev/null +++ b/Eigen/src/Sparse/SparseDot.h @@ -0,0 +1,97 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_SPARSE_DOT_H +#define EIGEN_SPARSE_DOT_H + +template<typename Derived> +template<typename OtherDerived> +typename ei_traits<Derived>::Scalar +SparseMatrixBase<Derived>::dot(const MatrixBase<OtherDerived>& other) const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived) + EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(Derived,OtherDerived) + EIGEN_STATIC_ASSERT((ei_is_same_type<Scalar, typename OtherDerived::Scalar>::ret), + YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) + + ei_assert(size() == other.size()); + ei_assert(other.size()>0 && "you are using a non initialized vector"); + + typename Derived::InnerIterator i(derived(),0); + Scalar res = 0; + while (i) + { + res += i.value() * ei_conj(other.coeff(i.index())); + ++i; + } + return res; +} + +template<typename Derived> +template<typename OtherDerived> +typename ei_traits<Derived>::Scalar +SparseMatrixBase<Derived>::dot(const SparseMatrixBase<OtherDerived>& other) const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived) + EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(Derived,OtherDerived) + EIGEN_STATIC_ASSERT((ei_is_same_type<Scalar, typename OtherDerived::Scalar>::ret), + YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) + + ei_assert(size() == other.size()); + + typename Derived::InnerIterator i(derived(),0); + typename OtherDerived::InnerIterator j(other.derived(),0); + Scalar res = 0; + while (i && j) + { + if (i.index()==j.index()) + { + res += i.value() * ei_conj(j.value()); + ++i; ++j; + } + else if (i.index()<j.index()) + ++i; + else + ++j; + } + return res; +} + +template<typename Derived> +inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real +SparseMatrixBase<Derived>::squaredNorm() const +{ + return ei_real((*this).cwise().abs2().sum()); +} + +template<typename Derived> +inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real +SparseMatrixBase<Derived>::norm() const +{ + return ei_sqrt(squaredNorm()); +} + +#endif // EIGEN_SPARSE_DOT_H diff --git a/Eigen/src/Sparse/SparseFlagged.h b/Eigen/src/Sparse/SparseFlagged.h new file mode 100644 index 000000000..2edac9c1e --- /dev/null +++ b/Eigen/src/Sparse/SparseFlagged.h @@ -0,0 +1,98 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_SPARSE_FLAGGED_H +#define EIGEN_SPARSE_FLAGGED_H + +template<typename ExpressionType, unsigned int Added, unsigned int Removed> +struct ei_traits<SparseFlagged<ExpressionType, Added, Removed> > : ei_traits<ExpressionType> +{ + enum { Flags = (ExpressionType::Flags | Added) & ~Removed }; +}; + +template<typename ExpressionType, unsigned int Added, unsigned int Removed> class SparseFlagged + : public SparseMatrixBase<SparseFlagged<ExpressionType, Added, Removed> > +{ + public: + + EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseFlagged) + class InnerIterator; + class ReverseInnerIterator; + + typedef typename ei_meta_if<ei_must_nest_by_value<ExpressionType>::ret, + ExpressionType, const ExpressionType&>::ret ExpressionTypeNested; + typedef typename ExpressionType::InnerIterator InnerIterator; + + inline SparseFlagged(const ExpressionType& matrix) : m_matrix(matrix) {} + + inline int rows() const { return m_matrix.rows(); } + inline int cols() const { return m_matrix.cols(); } + + // FIXME should be keep them ? + inline Scalar& coeffRef(int row, int col) + { return m_matrix.const_cast_derived().coeffRef(col, row); } + + inline const Scalar coeff(int row, int col) const + { return m_matrix.coeff(col, row); } + + inline const Scalar coeff(int index) const + { return m_matrix.coeff(index); } + + inline Scalar& coeffRef(int index) + { return m_matrix.const_cast_derived().coeffRef(index); } + + protected: + ExpressionTypeNested m_matrix; +}; + +template<typename ExpressionType, unsigned int Added, unsigned int Removed> + class SparseFlagged<ExpressionType,Added,Removed>::InnerIterator : public ExpressionType::InnerIterator +{ + public: + + EIGEN_STRONG_INLINE InnerIterator(const SparseFlagged& xpr, int outer) + : ExpressionType::InnerIterator(xpr.m_matrix, outer) + {} +}; + +template<typename ExpressionType, unsigned int Added, unsigned int Removed> + class SparseFlagged<ExpressionType,Added,Removed>::ReverseInnerIterator : public ExpressionType::ReverseInnerIterator +{ + public: + + EIGEN_STRONG_INLINE ReverseInnerIterator(const SparseFlagged& xpr, int outer) + : ExpressionType::ReverseInnerIterator(xpr.m_matrix, outer) + {} +}; + +template<typename Derived> +template<unsigned int Added> +inline const SparseFlagged<Derived, Added, 0> +SparseMatrixBase<Derived>::marked() const +{ + return derived(); +} + +#endif // EIGEN_SPARSE_FLAGGED_H diff --git a/Eigen/src/Sparse/SparseFuzzy.h b/Eigen/src/Sparse/SparseFuzzy.h new file mode 100644 index 000000000..355f4d52e --- /dev/null +++ b/Eigen/src/Sparse/SparseFuzzy.h @@ -0,0 +1,41 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_SPARSE_FUZZY_H +#define EIGEN_SPARSE_FUZZY_H + +// template<typename Derived> +// template<typename OtherDerived> +// bool SparseMatrixBase<Derived>::isApprox( +// const OtherDerived& other, +// typename NumTraits<Scalar>::Real prec +// ) const +// { +// const typename ei_nested<Derived,2>::type nested(derived()); +// const typename ei_nested<OtherDerived,2>::type otherNested(other.derived()); +// return (nested - otherNested).cwise().abs2().sum() +// <= prec * prec * std::min(nested.cwise().abs2().sum(), otherNested.cwise().abs2().sum()); +// } + +#endif // EIGEN_SPARSE_FUZZY_H diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h index 6820ae403..a732bdc31 100644 --- a/Eigen/src/Sparse/SparseMatrix.h +++ b/Eigen/src/Sparse/SparseMatrix.h @@ -56,14 +56,14 @@ class SparseMatrix : public SparseMatrixBase<SparseMatrix<_Scalar, _Flags> > { public: - EIGEN_GENERIC_PUBLIC_INTERFACE(SparseMatrix) + EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseMatrix) protected: public: typedef SparseMatrixBase<SparseMatrix> SparseBase; enum { - RowMajor = SparseBase::RowMajor + RowMajor = SparseBase::IsRowMajor }; typedef SparseMatrix<Scalar,(Flags&~RowMajorBit)|(RowMajor?RowMajorBit:0)> TransposedSparseMatrix; @@ -267,7 +267,7 @@ class SparseMatrix } template<typename OtherDerived> - inline SparseMatrix(const MatrixBase<OtherDerived>& other) + inline SparseMatrix(const SparseMatrixBase<OtherDerived>& other) : m_outerSize(0), m_innerSize(0), m_outerIndex(0) { *this = other.derived(); @@ -305,7 +305,7 @@ class SparseMatrix } template<typename OtherDerived> - inline SparseMatrix& operator=(const MatrixBase<OtherDerived>& other) + inline SparseMatrix& operator=(const SparseMatrixBase<OtherDerived>& other) { const bool needToTranspose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit); if (needToTranspose) diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h index d70c99259..d01fa1ec5 100644 --- a/Eigen/src/Sparse/SparseMatrixBase.h +++ b/Eigen/src/Sparse/SparseMatrixBase.h @@ -25,43 +25,131 @@ #ifndef EIGEN_SPARSEMATRIXBASE_H #define EIGEN_SPARSEMATRIXBASE_H -template<typename Derived> -class SparseMatrixBase : public MatrixBase<Derived> +template<typename Derived> class SparseMatrixBase { public: - typedef MatrixBase<Derived> Base; - typedef typename Base::Scalar Scalar; - typedef typename Base::RealScalar RealScalar; + typedef typename ei_traits<Derived>::Scalar Scalar; +// typedef typename Derived::InnerIterator InnerIterator; + enum { - Flags = Base::Flags, - RowMajor = ei_traits<Derived>::Flags&RowMajorBit ? 1 : 0 + + RowsAtCompileTime = ei_traits<Derived>::RowsAtCompileTime, + /**< The number of rows at compile-time. This is just a copy of the value provided + * by the \a Derived type. If a value is not known at compile-time, + * it is set to the \a Dynamic constant. + * \sa MatrixBase::rows(), MatrixBase::cols(), ColsAtCompileTime, SizeAtCompileTime */ + + ColsAtCompileTime = ei_traits<Derived>::ColsAtCompileTime, + /**< The number of columns at compile-time. This is just a copy of the value provided + * by the \a Derived type. If a value is not known at compile-time, + * it is set to the \a Dynamic constant. + * \sa MatrixBase::rows(), MatrixBase::cols(), RowsAtCompileTime, SizeAtCompileTime */ + + + SizeAtCompileTime = (ei_size_at_compile_time<ei_traits<Derived>::RowsAtCompileTime, + ei_traits<Derived>::ColsAtCompileTime>::ret), + /**< This is equal to the number of coefficients, i.e. the number of + * rows times the number of columns, or to \a Dynamic if this is not + * known at compile-time. \sa RowsAtCompileTime, ColsAtCompileTime */ + + MaxRowsAtCompileTime = RowsAtCompileTime, + MaxColsAtCompileTime = ColsAtCompileTime, + + MaxSizeAtCompileTime = (ei_size_at_compile_time<MaxRowsAtCompileTime, + MaxColsAtCompileTime>::ret), + + IsVectorAtCompileTime = RowsAtCompileTime == 1 || ColsAtCompileTime == 1, + /**< This is set to true if either the number of rows or the number of + * columns is known at compile-time to be equal to 1. Indeed, in that case, + * we are dealing with a column-vector (if there is only one column) or with + * a row-vector (if there is only one row). */ + + Flags = ei_traits<Derived>::Flags, + /**< This stores expression \ref flags flags which may or may not be inherited by new expressions + * constructed from this one. See the \ref flags "list of flags". + */ + + CoeffReadCost = ei_traits<Derived>::CoeffReadCost, + /**< This is a rough measure of how expensive it is to read one coefficient from + * this expression. + */ + + IsRowMajor = Flags&RowMajorBit ? 1 : 0 }; + /** \internal the return type of MatrixBase::conjugate() */ + typedef typename ei_meta_if<NumTraits<Scalar>::IsComplex, + const SparseCwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, Derived>, + const Derived& + >::ret ConjugateReturnType; + /** \internal the return type of MatrixBase::real() */ + typedef CwiseUnaryOp<ei_scalar_real_op<Scalar>, Derived> RealReturnType; + /** \internal the return type of MatrixBase::imag() */ + typedef CwiseUnaryOp<ei_scalar_imag_op<Scalar>, Derived> ImagReturnType; + /** \internal the return type of MatrixBase::adjoint() */ + typedef Eigen::Transpose<NestByValue<typename ei_cleantype<ConjugateReturnType>::type> > + AdjointReturnType; + +#ifndef EIGEN_PARSED_BY_DOXYGEN + /** This is the "real scalar" type; if the \a Scalar type is already real numbers + * (e.g. int, float or double) then \a RealScalar is just the same as \a Scalar. If + * \a Scalar is \a std::complex<T> then RealScalar is \a T. + * + * \sa class NumTraits + */ + typedef typename NumTraits<Scalar>::Real RealScalar; + + /** type of the equivalent square matrix */ + typedef Matrix<Scalar,EIGEN_ENUM_MAX(RowsAtCompileTime,ColsAtCompileTime), + EIGEN_ENUM_MAX(RowsAtCompileTime,ColsAtCompileTime)> SquareMatrixType; + inline const Derived& derived() const { return *static_cast<const Derived*>(this); } inline Derived& derived() { return *static_cast<Derived*>(this); } inline Derived& const_cast_derived() const { return *static_cast<Derived*>(const_cast<SparseMatrixBase*>(this)); } +#endif // not EIGEN_PARSED_BY_DOXYGEN - SparseMatrixBase() - : m_isRValue(false) - {} + /** \returns the number of rows. \sa cols(), RowsAtCompileTime */ + inline int rows() const { return derived().rows(); } + /** \returns the number of columns. \sa rows(), ColsAtCompileTime*/ + inline int cols() const { return derived().cols(); } + /** \returns the number of coefficients, which is \a rows()*cols(). + * \sa rows(), cols(), SizeAtCompileTime. */ + inline int size() const { return rows() * cols(); } + /** \returns the number of nonzero coefficients which is in practice the number + * of stored coefficients. */ + inline int nonZeros() const { return derived.nonZeros(); } + /** \returns true if either the number of rows or the number of columns is equal to 1. + * In other words, this function returns + * \code rows()==1 || cols()==1 \endcode + * \sa rows(), cols(), IsVectorAtCompileTime. */ + inline bool isVector() const { return rows()==1 || cols()==1; } + /** \returns the size of the storage major dimension, + * i.e., the number of columns for a columns major matrix, and the number of rows otherwise */ + int outerSize() const { return (int(Flags)&RowMajorBit) ? this->rows() : this->cols(); } + /** \returns the size of the inner dimension according to the storage order, + * i.e., the number of rows for a columns major matrix, and the number of cols otherwise */ + int innerSize() const { return (int(Flags)&RowMajorBit) ? this->cols() : this->rows(); } bool isRValue() const { return m_isRValue; } Derived& markAsRValue() { m_isRValue = true; return derived(); } + SparseMatrixBase() : m_isRValue(false) { /* TODO check flags */ } + inline Derived& operator=(const Derived& other) { // std::cout << "Derived& operator=(const Derived& other)\n"; - if (other.isRValue()) - derived().swap(other.const_cast_derived()); - else +// if (other.isRValue()) +// derived().swap(other.const_cast_derived()); +// else this->operator=<Derived>(other); return derived(); } + template<typename OtherDerived> - inline Derived& operator=(const MatrixBase<OtherDerived>& other) + inline void assignGeneric(const OtherDerived& other) { // std::cout << "Derived& operator=(const MatrixBase<OtherDerived>& other)\n"; //const bool transpose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit); @@ -87,9 +175,9 @@ class SparseMatrixBase : public MatrixBase<Derived> temp.endFill(); derived() = temp.markAsRValue(); - return derived(); } + template<typename OtherDerived> inline Derived& operator=(const SparseMatrixBase<OtherDerived>& other) { @@ -110,7 +198,7 @@ class SparseMatrixBase : public MatrixBase<Derived> Scalar v = it.value(); if (v!=Scalar(0)) { - if (RowMajor) derived().fill(j,it.index()) = v; + if (IsRowMajor) derived().fill(j,it.index()) = v; else derived().fill(it.index(),j) = v; } } @@ -118,12 +206,14 @@ class SparseMatrixBase : public MatrixBase<Derived> derived().endFill(); } else - this->operator=<OtherDerived>(static_cast<const MatrixBase<OtherDerived>&>(other)); + { + assignGeneric(other.derived()); + } return derived(); } template<typename Lhs, typename Rhs> - inline Derived& operator=(const Product<Lhs,Rhs,SparseProduct>& product); + inline Derived& operator=(const SparseProduct<Lhs,Rhs>& product); friend std::ostream & operator << (std::ostream & s, const SparseMatrixBase& m) { @@ -167,6 +257,292 @@ class SparseMatrixBase : public MatrixBase<Derived> return s; } + const SparseCwiseUnaryOp<ei_scalar_opposite_op<typename ei_traits<Derived>::Scalar>,Derived> operator-() const; + + template<typename OtherDerived> + const SparseCwiseBinaryOp<ei_scalar_sum_op<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived> + operator+(const SparseMatrixBase<OtherDerived> &other) const; + + template<typename OtherDerived> + const SparseCwiseBinaryOp<ei_scalar_difference_op<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived> + operator-(const SparseMatrixBase<OtherDerived> &other) const; + + template<typename OtherDerived> + Derived& operator+=(const SparseMatrixBase<OtherDerived>& other); + template<typename OtherDerived> + Derived& operator-=(const SparseMatrixBase<OtherDerived>& other); + +// template<typename Lhs,typename Rhs> +// Derived& operator+=(const Flagged<Product<Lhs,Rhs,CacheFriendlyProduct>, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit>& other); + + Derived& operator*=(const Scalar& other); + Derived& operator/=(const Scalar& other); + + const SparseCwiseUnaryOp<ei_scalar_multiple_op<typename ei_traits<Derived>::Scalar>, Derived> + operator*(const Scalar& scalar) const; + const SparseCwiseUnaryOp<ei_scalar_quotient1_op<typename ei_traits<Derived>::Scalar>, Derived> + operator/(const Scalar& scalar) const; + + inline friend const SparseCwiseUnaryOp<ei_scalar_multiple_op<typename ei_traits<Derived>::Scalar>, Derived> + operator*(const Scalar& scalar, const SparseMatrixBase& matrix) + { return matrix*scalar; } + + + template<typename OtherDerived> + const typename SparseProductReturnType<Derived,OtherDerived>::Type + operator*(const SparseMatrixBase<OtherDerived> &other) const; + + template<typename OtherDerived> + Derived& operator*=(const SparseMatrixBase<OtherDerived>& other); + + template<typename OtherDerived> + typename ei_plain_matrix_type_column_major<OtherDerived>::type + solveTriangular(const MatrixBase<OtherDerived>& other) const; + + template<typename OtherDerived> + void solveTriangularInPlace(MatrixBase<OtherDerived>& other) const; + + template<typename OtherDerived> Scalar dot(const MatrixBase<OtherDerived>& other) const; + template<typename OtherDerived> Scalar dot(const SparseMatrixBase<OtherDerived>& other) const; + RealScalar squaredNorm() const; + RealScalar norm() const; +// const PlainMatrixType normalized() const; +// void normalize(); + + SparseTranspose<Derived> transpose() { return derived(); } + const SparseTranspose<Derived> transpose() const { return derived(); } + // void transposeInPlace(); + // const AdjointReturnType adjoint() const; + + SparseInnerVector<Derived> innerVector(int outer); + const SparseInnerVector<Derived> innerVector(int outer) const; + +// RowXpr row(int i); +// const RowXpr row(int i) const; + +// ColXpr col(int i); +// const ColXpr col(int i) const; + +// typename BlockReturnType<Derived>::Type block(int startRow, int startCol, int blockRows, int blockCols); +// const typename BlockReturnType<Derived>::Type +// block(int startRow, int startCol, int blockRows, int blockCols) const; +// +// typename BlockReturnType<Derived>::SubVectorType segment(int start, int size); +// const typename BlockReturnType<Derived>::SubVectorType segment(int start, int size) const; +// +// typename BlockReturnType<Derived,Dynamic>::SubVectorType start(int size); +// const typename BlockReturnType<Derived,Dynamic>::SubVectorType start(int size) const; +// +// typename BlockReturnType<Derived,Dynamic>::SubVectorType end(int size); +// const typename BlockReturnType<Derived,Dynamic>::SubVectorType end(int size) const; +// +// typename BlockReturnType<Derived>::Type corner(CornerType type, int cRows, int cCols); +// const typename BlockReturnType<Derived>::Type corner(CornerType type, int cRows, int cCols) const; +// +// template<int BlockRows, int BlockCols> +// typename BlockReturnType<Derived, BlockRows, BlockCols>::Type block(int startRow, int startCol); +// template<int BlockRows, int BlockCols> +// const typename BlockReturnType<Derived, BlockRows, BlockCols>::Type block(int startRow, int startCol) const; + +// template<int CRows, int CCols> +// typename BlockReturnType<Derived, CRows, CCols>::Type corner(CornerType type); +// template<int CRows, int CCols> +// const typename BlockReturnType<Derived, CRows, CCols>::Type corner(CornerType type) const; + +// template<int Size> typename BlockReturnType<Derived,Size>::SubVectorType start(void); +// template<int Size> const typename BlockReturnType<Derived,Size>::SubVectorType start() const; + +// template<int Size> typename BlockReturnType<Derived,Size>::SubVectorType end(); +// template<int Size> const typename BlockReturnType<Derived,Size>::SubVectorType end() const; + +// template<int Size> typename BlockReturnType<Derived,Size>::SubVectorType segment(int start); +// template<int Size> const typename BlockReturnType<Derived,Size>::SubVectorType segment(int start) const; + +// DiagonalCoeffs<Derived> diagonal(); +// const DiagonalCoeffs<Derived> diagonal() const; + +// template<unsigned int Mode> Part<Derived, Mode> part(); +// template<unsigned int Mode> const Part<Derived, Mode> part() const; + + +// static const ConstantReturnType Constant(int rows, int cols, const Scalar& value); +// static const ConstantReturnType Constant(int size, const Scalar& value); +// static const ConstantReturnType Constant(const Scalar& value); + +// template<typename CustomNullaryOp> +// static const CwiseNullaryOp<CustomNullaryOp, Derived> NullaryExpr(int rows, int cols, const CustomNullaryOp& func); +// template<typename CustomNullaryOp> +// static const CwiseNullaryOp<CustomNullaryOp, Derived> NullaryExpr(int size, const CustomNullaryOp& func); +// template<typename CustomNullaryOp> +// static const CwiseNullaryOp<CustomNullaryOp, Derived> NullaryExpr(const CustomNullaryOp& func); + +// static const ConstantReturnType Zero(int rows, int cols); +// static const ConstantReturnType Zero(int size); +// static const ConstantReturnType Zero(); +// static const ConstantReturnType Ones(int rows, int cols); +// static const ConstantReturnType Ones(int size); +// static const ConstantReturnType Ones(); +// static const IdentityReturnType Identity(); +// static const IdentityReturnType Identity(int rows, int cols); +// static const BasisReturnType Unit(int size, int i); +// static const BasisReturnType Unit(int i); +// static const BasisReturnType UnitX(); +// static const BasisReturnType UnitY(); +// static const BasisReturnType UnitZ(); +// static const BasisReturnType UnitW(); + +// const DiagonalMatrix<Derived> asDiagonal() const; + +// Derived& setConstant(const Scalar& value); +// Derived& setZero(); +// Derived& setOnes(); +// Derived& setRandom(); +// Derived& setIdentity(); + + Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime> toDense() const + { + Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime> res(rows(),cols()); + res.setZero(); + for (int j=0; j<outerSize(); ++j) + { + for (typename Derived::InnerIterator i(derived(),j); i; ++i) + if(IsRowMajor) + res.coeffRef(j,i.index()) = i.value(); + else + res.coeffRef(i.index(),j) = i.value(); + } + return res; + } + + template<typename OtherDerived> + bool isApprox(const SparseMatrixBase<OtherDerived>& other, + RealScalar prec = precision<Scalar>()) const + { return toDense().isApprox(other.toDense(),prec); } + + template<typename OtherDerived> + bool isApprox(const MatrixBase<OtherDerived>& other, + RealScalar prec = precision<Scalar>()) const + { return toDense().isApprox(other,prec); } +// bool isMuchSmallerThan(const RealScalar& other, +// RealScalar prec = precision<Scalar>()) const; +// template<typename OtherDerived> +// bool isMuchSmallerThan(const MatrixBase<OtherDerived>& other, +// RealScalar prec = precision<Scalar>()) const; + +// bool isApproxToConstant(const Scalar& value, RealScalar prec = precision<Scalar>()) const; +// bool isZero(RealScalar prec = precision<Scalar>()) const; +// bool isOnes(RealScalar prec = precision<Scalar>()) const; +// bool isIdentity(RealScalar prec = precision<Scalar>()) const; +// bool isDiagonal(RealScalar prec = precision<Scalar>()) const; + +// bool isUpperTriangular(RealScalar prec = precision<Scalar>()) const; +// bool isLowerTriangular(RealScalar prec = precision<Scalar>()) const; + +// template<typename OtherDerived> +// bool isOrthogonal(const MatrixBase<OtherDerived>& other, +// RealScalar prec = precision<Scalar>()) const; +// bool isUnitary(RealScalar prec = precision<Scalar>()) const; + +// template<typename OtherDerived> +// inline bool operator==(const MatrixBase<OtherDerived>& other) const +// { return (cwise() == other).all(); } + +// template<typename OtherDerived> +// inline bool operator!=(const MatrixBase<OtherDerived>& other) const +// { return (cwise() != other).any(); } + + + template<typename NewType> + const SparseCwiseUnaryOp<ei_scalar_cast_op<typename ei_traits<Derived>::Scalar, NewType>, Derived> cast() const; + + /** \returns the matrix or vector obtained by evaluating this expression. + * + * Notice that in the case of a plain matrix or vector (not an expression) this function just returns + * a const reference, in order to avoid a useless copy. + */ + EIGEN_STRONG_INLINE const typename ei_eval<Derived>::type eval() const + { return typename ei_eval<Derived>::type(derived()); } + +// template<typename OtherDerived> +// void swap(const MatrixBase<OtherDerived>& other); + + template<unsigned int Added> + const SparseFlagged<Derived, Added, 0> marked() const; +// const Flagged<Derived, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit> lazy() const; + + /** \returns number of elements to skip to pass from one row (resp. column) to another + * for a row-major (resp. column-major) matrix. + * Combined with coeffRef() and the \ref flags flags, it allows a direct access to the data + * of the underlying matrix. + */ +// inline int stride(void) const { return derived().stride(); } + +// inline const NestByValue<Derived> nestByValue() const; + + + ConjugateReturnType conjugate() const; + const RealReturnType real() const; + const ImagReturnType imag() const; + + template<typename CustomUnaryOp> + const SparseCwiseUnaryOp<CustomUnaryOp, Derived> unaryExpr(const CustomUnaryOp& func = CustomUnaryOp()) const; + +// template<typename CustomBinaryOp, typename OtherDerived> +// const CwiseBinaryOp<CustomBinaryOp, Derived, OtherDerived> +// binaryExpr(const MatrixBase<OtherDerived> &other, const CustomBinaryOp& func = CustomBinaryOp()) const; + + + Scalar sum() const; +// Scalar trace() const; + +// typename ei_traits<Derived>::Scalar minCoeff() const; +// typename ei_traits<Derived>::Scalar maxCoeff() const; + +// typename ei_traits<Derived>::Scalar minCoeff(int* row, int* col = 0) const; +// typename ei_traits<Derived>::Scalar maxCoeff(int* row, int* col = 0) const; + +// template<typename BinaryOp> +// typename ei_result_of<BinaryOp(typename ei_traits<Derived>::Scalar)>::type +// redux(const BinaryOp& func) const; + +// template<typename Visitor> +// void visit(Visitor& func) const; + + + const SparseCwise<Derived> cwise() const; + SparseCwise<Derived> cwise(); + +// inline const WithFormat<Derived> format(const IOFormat& fmt) const; + +/////////// Array module /////////// + /* + bool all(void) const; + bool any(void) const; + + const PartialRedux<Derived,Horizontal> rowwise() const; + const PartialRedux<Derived,Vertical> colwise() const; + + static const CwiseNullaryOp<ei_scalar_random_op<Scalar>,Derived> Random(int rows, int cols); + static const CwiseNullaryOp<ei_scalar_random_op<Scalar>,Derived> Random(int size); + static const CwiseNullaryOp<ei_scalar_random_op<Scalar>,Derived> Random(); + + template<typename ThenDerived,typename ElseDerived> + const Select<Derived,ThenDerived,ElseDerived> + select(const MatrixBase<ThenDerived>& thenMatrix, + const MatrixBase<ElseDerived>& elseMatrix) const; + + template<typename ThenDerived> + inline const Select<Derived,ThenDerived, NestByValue<typename ThenDerived::ConstantReturnType> > + select(const MatrixBase<ThenDerived>& thenMatrix, typename ThenDerived::Scalar elseScalar) const; + + template<typename ElseDerived> + inline const Select<Derived, NestByValue<typename ElseDerived::ConstantReturnType>, ElseDerived > + select(typename ElseDerived::Scalar thenScalar, const MatrixBase<ElseDerived>& elseMatrix) const; + + template<int p> RealScalar lpNorm() const; + */ + + // template<typename OtherDerived> // Scalar dot(const MatrixBase<OtherDerived>& other) const // { diff --git a/Eigen/src/Sparse/SparseProduct.h b/Eigen/src/Sparse/SparseProduct.h index 36dfb31ac..b4ba2ee6f 100644 --- a/Eigen/src/Sparse/SparseProduct.h +++ b/Eigen/src/Sparse/SparseProduct.h @@ -27,7 +27,7 @@ // sparse product return type specialization template<typename Lhs, typename Rhs> -struct ProductReturnType<Lhs,Rhs,SparseProduct> +struct SparseProductReturnType { typedef typename ei_traits<Lhs>::Scalar Scalar; enum { @@ -47,12 +47,11 @@ struct ProductReturnType<Lhs,Rhs,SparseProduct> SparseMatrix<Scalar,0>, const typename ei_nested<Rhs,Lhs::RowsAtCompileTime>::type>::ret RhsNested; - typedef Product<LhsNested, - RhsNested, SparseProduct> Type; + typedef SparseProduct<LhsNested, RhsNested> Type; }; template<typename LhsNested, typename RhsNested> -struct ei_traits<Product<LhsNested, RhsNested, SparseProduct> > +struct ei_traits<SparseProduct<LhsNested, RhsNested> > { // clean the nested types: typedef typename ei_cleantype<LhsNested>::type _LhsNested; @@ -87,30 +86,28 @@ struct ei_traits<Product<LhsNested, RhsNested, SparseProduct> > }; }; -template<typename LhsNested, typename RhsNested> class Product<LhsNested,RhsNested,SparseProduct> : ei_no_assignment_operator, - public MatrixBase<Product<LhsNested, RhsNested, SparseProduct> > +template<typename LhsNested, typename RhsNested> +class SparseProduct : ei_no_assignment_operator, + public SparseMatrixBase<SparseProduct<LhsNested, RhsNested> > { public: - EIGEN_GENERIC_PUBLIC_INTERFACE(Product) + EIGEN_GENERIC_PUBLIC_INTERFACE(SparseProduct) private: - typedef typename ei_traits<Product>::_LhsNested _LhsNested; - typedef typename ei_traits<Product>::_RhsNested _RhsNested; + typedef typename ei_traits<SparseProduct>::_LhsNested _LhsNested; + typedef typename ei_traits<SparseProduct>::_RhsNested _RhsNested; public: template<typename Lhs, typename Rhs> - inline Product(const Lhs& lhs, const Rhs& rhs) + inline SparseProduct(const Lhs& lhs, const Rhs& rhs) : m_lhs(lhs), m_rhs(rhs) { ei_assert(lhs.cols() == rhs.rows()); } - Scalar coeff(int, int) const { ei_assert(false && "eigen internal error"); } - Scalar& coeffRef(int, int) { ei_assert(false && "eigen internal error"); } - inline int rows() const { return m_lhs.rows(); } inline int cols() const { return m_rhs.cols(); } @@ -231,28 +228,51 @@ struct ei_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,ColMajor> // by ProductReturnType which transform it to (col col *) by evaluating rhs. +// template<typename Derived> +// template<typename Lhs, typename Rhs> +// inline Derived& SparseMatrixBase<Derived>::lazyAssign(const SparseProduct<Lhs,Rhs>& product) +// { +// // std::cout << "sparse product to dense\n"; +// ei_sparse_product_selector< +// typename ei_cleantype<Lhs>::type, +// typename ei_cleantype<Rhs>::type, +// typename ei_cleantype<Derived>::type>::run(product.lhs(),product.rhs(),derived()); +// return derived(); +// } + template<typename Derived> template<typename Lhs, typename Rhs> -inline Derived& MatrixBase<Derived>::lazyAssign(const Product<Lhs,Rhs,SparseProduct>& product) +inline Derived& SparseMatrixBase<Derived>::operator=(const SparseProduct<Lhs,Rhs>& product) { -// std::cout << "sparse product to dense\n"; +// std::cout << "sparse product to sparse\n"; ei_sparse_product_selector< typename ei_cleantype<Lhs>::type, typename ei_cleantype<Rhs>::type, - typename ei_cleantype<Derived>::type>::run(product.lhs(),product.rhs(),derived()); + Derived>::run(product.lhs(),product.rhs(),derived()); return derived(); } template<typename Derived> -template<typename Lhs, typename Rhs> -inline Derived& SparseMatrixBase<Derived>::operator=(const Product<Lhs,Rhs,SparseProduct>& product) +template<typename OtherDerived> +inline const typename SparseProductReturnType<Derived,OtherDerived>::Type +SparseMatrixBase<Derived>::operator*(const SparseMatrixBase<OtherDerived> &other) const { -// std::cout << "sparse product to sparse\n"; - ei_sparse_product_selector< - typename ei_cleantype<Lhs>::type, - typename ei_cleantype<Rhs>::type, - Derived>::run(product.lhs(),product.rhs(),derived()); - return derived(); + enum { + ProductIsValid = Derived::ColsAtCompileTime==Dynamic + || OtherDerived::RowsAtCompileTime==Dynamic + || int(Derived::ColsAtCompileTime)==int(OtherDerived::RowsAtCompileTime), + AreVectors = Derived::IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime, + SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(Derived,OtherDerived) + }; + // note to the lost user: + // * for a dot product use: v1.dot(v2) + // * for a coeff-wise product use: v1.cwise()*v2 + EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes), + INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS) + EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors), + INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION) + EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT) + return typename SparseProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived()); } #endif // EIGEN_SPARSEPRODUCT_H diff --git a/Eigen/src/Sparse/SparseRedux.h b/Eigen/src/Sparse/SparseRedux.h index 746ae2f8a..f0d370548 100644 --- a/Eigen/src/Sparse/SparseRedux.h +++ b/Eigen/src/Sparse/SparseRedux.h @@ -25,46 +25,16 @@ #ifndef EIGEN_SPARSEREDUX_H #define EIGEN_SPARSEREDUX_H - -template<typename Derived, int Vectorization, int Unrolling> -struct ei_sum_impl<Derived, Vectorization, Unrolling, IsSparse> -{ - typedef typename Derived::Scalar Scalar; - static Scalar run(const Derived& mat) - { - ei_assert(mat.rows()>0 && mat.cols()>0 && "you are using a non initialized matrix"); - Scalar res = 0; - for (int j=0; j<mat.outerSize(); ++j) - for (typename Derived::InnerIterator iter(mat,j); iter; ++iter) - res += iter.value(); - return res; - } -}; - -template<typename Derived1, typename Derived2, int Vectorization, int Unrolling> -struct ei_dot_impl<Derived1, Derived2, Vectorization, Unrolling, IsSparse> +template<typename Derived> +typename ei_traits<Derived>::Scalar +SparseMatrixBase<Derived>::sum() const { - typedef typename Derived1::Scalar Scalar; - static Scalar run(const Derived1& v1, const Derived2& v2) - { - ei_assert(v1.size()>0 && "you are using a non initialized vector"); - typename Derived1::InnerIterator i(v1,0); - typename Derived2::InnerIterator j(v2,0); - Scalar res = 0; - while (i && j) - { - if (i.index()==j.index()) - { - res += i.value() * ei_conj(j.value()); - ++i; ++j; - } - else if (i.index()<j.index()) - ++i; - else - ++j; - } - return res; - } -}; + ei_assert(rows()>0 && cols()>0 && "you are using a non initialized matrix"); + Scalar res = 0; + for (int j=0; j<outerSize(); ++j) + for (typename Derived::InnerIterator iter(derived(),j); iter; ++iter) + res += iter.value(); + return res; +} #endif // EIGEN_SPARSEREDUX_H diff --git a/Eigen/src/Sparse/SparseTranspose.h b/Eigen/src/Sparse/SparseTranspose.h new file mode 100644 index 000000000..89a14d707 --- /dev/null +++ b/Eigen/src/Sparse/SparseTranspose.h @@ -0,0 +1,85 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_SPARSETRANSPOSE_H +#define EIGEN_SPARSETRANSPOSE_H + +template<typename MatrixType> +struct ei_traits<SparseTranspose<MatrixType> > : ei_traits<Transpose<MatrixType> > +{}; + +template<typename MatrixType> class SparseTranspose + : public SparseMatrixBase<SparseTranspose<MatrixType> > +{ + public: + + EIGEN_GENERIC_PUBLIC_INTERFACE(SparseTranspose) + + class InnerIterator; + class ReverseInnerIterator; + + inline SparseTranspose(const MatrixType& matrix) : m_matrix(matrix) {} + + //EIGEN_INHERIT_ASSIGNMENT_OPERATORS(SparseTranspose) + + inline int rows() const { return m_matrix.cols(); } + inline int cols() const { return m_matrix.rows(); } + inline int nonZeros() const { return m_matrix.nonZeros(); } + + // FIXME should be keep them ? + inline Scalar& coeffRef(int row, int col) + { return m_matrix.const_cast_derived().coeffRef(col, row); } + + inline const Scalar coeff(int row, int col) const + { return m_matrix.coeff(col, row); } + + inline const Scalar coeff(int index) const + { return m_matrix.coeff(index); } + + inline Scalar& coeffRef(int index) + { return m_matrix.const_cast_derived().coeffRef(index); } + + protected: + const typename MatrixType::Nested m_matrix; +}; + +template<typename MatrixType> class SparseTranspose<MatrixType>::InnerIterator : public MatrixType::InnerIterator +{ + public: + + EIGEN_STRONG_INLINE InnerIterator(const SparseTranspose& trans, int outer) + : MatrixType::InnerIterator(trans.m_matrix, outer) + {} +}; + +template<typename MatrixType> class SparseTranspose<MatrixType>::ReverseInnerIterator : public MatrixType::ReverseInnerIterator +{ + public: + + EIGEN_STRONG_INLINE ReverseInnerIterator(const SparseTranspose& xpr, int outer) + : MatrixType::ReverseInnerIterator(xpr.m_matrix, outer) + {} +}; + +#endif // EIGEN_SPARSETRANSPOSE_H diff --git a/Eigen/src/Sparse/SparseUtil.h b/Eigen/src/Sparse/SparseUtil.h index dcafcc008..724fb9efb 100644 --- a/Eigen/src/Sparse/SparseUtil.h +++ b/Eigen/src/Sparse/SparseUtil.h @@ -31,6 +31,46 @@ #define EIGEN_DBG_SPARSE(X) X #endif +#define EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(Derived, Op) \ +template<typename OtherDerived> \ +EIGEN_STRONG_INLINE Derived& operator Op(const Eigen::SparseMatrixBase<OtherDerived>& other) \ +{ \ + return Base::operator Op(other.derived()); \ +} \ +EIGEN_STRONG_INLINE Derived& operator Op(const Derived& other) \ +{ \ + return Base::operator Op(other); \ +} + +#define EIGEN_SPARSE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, Op) \ +template<typename Other> \ +EIGEN_STRONG_INLINE Derived& operator Op(const Other& scalar) \ +{ \ + return Base::operator Op(scalar); \ +} + +#define EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATORS(Derived) \ +EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(Derived, =) \ +EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(Derived, +=) \ +EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(Derived, -=) \ +EIGEN_SPARSE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, *=) \ +EIGEN_SPARSE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, /=) + +#define _EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(Derived, BaseClass) \ +typedef BaseClass Base; \ +typedef typename Eigen::ei_traits<Derived>::Scalar Scalar; \ +typedef typename Eigen::NumTraits<Scalar>::Real RealScalar; \ +typedef typename Eigen::ei_nested<Derived>::type Nested; \ +enum { RowsAtCompileTime = Eigen::ei_traits<Derived>::RowsAtCompileTime, \ + ColsAtCompileTime = Eigen::ei_traits<Derived>::ColsAtCompileTime, \ + Flags = Eigen::ei_traits<Derived>::Flags, \ + CoeffReadCost = Eigen::ei_traits<Derived>::CoeffReadCost, \ + SizeAtCompileTime = Base::SizeAtCompileTime, \ + IsVectorAtCompileTime = Base::IsVectorAtCompileTime }; + +#define EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(Derived) \ +_EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(Derived, Eigen::SparseMatrixBase<Derived>) + enum SparseBackend { DefaultBackend, Taucs, @@ -64,6 +104,16 @@ template<typename Derived> class SparseMatrixBase; template<typename _Scalar, int _Flags = 0> class SparseMatrix; template<typename _Scalar, int _Flags = 0> class SparseVector; +template<typename MatrixType> class SparseTranspose; +template<typename MatrixType> class SparseInnerVector; +template<typename Derived> class SparseCwise; +template<typename UnaryOp, typename MatrixType> class SparseCwiseUnaryOp; +template<typename BinaryOp, typename Lhs, typename Rhs> class SparseCwiseBinaryOp; +template<typename Lhs, typename Rhs> class SparseProduct; +template<typename ExpressionType, unsigned int Added, unsigned int Removed> class SparseFlagged; + +template<typename Lhs, typename Rhs> struct SparseProductReturnType; + const int AccessPatternNotSupported = 0x0; const int AccessPatternSupported = 0x1; diff --git a/Eigen/src/Sparse/SparseVector.h b/Eigen/src/Sparse/SparseVector.h index 7a2ce70e5..b189cb5ee 100644 --- a/Eigen/src/Sparse/SparseVector.h +++ b/Eigen/src/Sparse/SparseVector.h @@ -58,7 +58,7 @@ class SparseVector : public SparseMatrixBase<SparseVector<_Scalar, _Flags> > { public: - EIGEN_GENERIC_PUBLIC_INTERFACE(SparseVector) + EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseVector) protected: public: diff --git a/Eigen/src/Sparse/TriangularSolver.h b/Eigen/src/Sparse/TriangularSolver.h index d44938273..8948ae45e 100644 --- a/Eigen/src/Sparse/TriangularSolver.h +++ b/Eigen/src/Sparse/TriangularSolver.h @@ -144,4 +144,35 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpperTriangular,ColMajor|IsSparse> } }; +template<typename Derived> +template<typename OtherDerived> +void SparseMatrixBase<Derived>::solveTriangularInPlace(MatrixBase<OtherDerived>& other) const +{ + ei_assert(derived().cols() == derived().rows()); + ei_assert(derived().cols() == other.rows()); + ei_assert(!(Flags & ZeroDiagBit)); + ei_assert(Flags & (UpperTriangularBit|LowerTriangularBit)); + + enum { copy = ei_traits<OtherDerived>::Flags & RowMajorBit }; + + typedef typename ei_meta_if<copy, + typename ei_plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::ret OtherCopy; + OtherCopy otherCopy(other.derived()); + + ei_solve_triangular_selector<Derived, typename ei_unref<OtherCopy>::type>::run(derived(), otherCopy); + + if (copy) + other = otherCopy; +} + +template<typename Derived> +template<typename OtherDerived> +typename ei_plain_matrix_type_column_major<OtherDerived>::type +SparseMatrixBase<Derived>::solveTriangular(const MatrixBase<OtherDerived>& other) const +{ + typename ei_plain_matrix_type_column_major<OtherDerived>::type res(other); + solveTriangularInPlace(res); + return res; +} + #endif // EIGEN_SPARSETRIANGULARSOLVER_H diff --git a/cmake/FindSuperLU.cmake b/cmake/FindSuperLU.cmake index e10561b60..fbefc943a 100644 --- a/cmake/FindSuperLU.cmake +++ b/cmake/FindSuperLU.cmake @@ -16,10 +16,6 @@ if(BLAS_FOUND) ) find_library(SUPERLU_LIBRARIES superlu PATHS $ENV{SUPERLUDIR} ${LIB_INSTALL_DIR}) - - if(SUPERLU_LIBRARIES AND CMAKE_COMPILER_IS_GNUCXX) - set(SUPERLU_LIBRARIES ${SUPERLU_LIBRARIES} -lgfortran) - endif(SUPERLU_LIBRARIES AND CMAKE_COMPILER_IS_GNUCXX) if(SUPERLU_LIBRARIES) set(SUPERLU_LIBRARIES ${SUPERLU_LIBRARIES} ${BLAS_LIBRARIES}) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index bc9a5ca73..021799f5e 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -195,8 +195,8 @@ ei_add_test(parametrizedline) ei_add_test(alignedbox) ei_add_test(regression) ei_add_test(stdvector) -ei_add_test(sparse_basic) ei_add_test(sparse_vector) +ei_add_test(sparse_basic) ei_add_test(sparse_solvers " " "${SPARSE_LIBS}") # print a summary of the different options diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp index c50682810..54272d871 100644 --- a/test/sparse_basic.cpp +++ b/test/sparse_basic.cpp @@ -72,7 +72,7 @@ template<typename Scalar> void sparse_basic(int rows, int cols) refMat.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5); VERIFY_IS_APPROX(m, refMat); -/* + /* // test InnerIterators and Block expressions for (int t=0; t<10; ++t) { @@ -81,23 +81,23 @@ template<typename Scalar> void sparse_basic(int rows, int cols) int w = ei_random<int>(1,cols-j-1); int h = ei_random<int>(1,rows-i-1); - VERIFY_IS_APPROX(m.block(i,j,h,w), refMat.block(i,j,h,w)); +// VERIFY_IS_APPROX(m.block(i,j,h,w), refMat.block(i,j,h,w)); for(int c=0; c<w; c++) { VERIFY_IS_APPROX(m.block(i,j,h,w).col(c), refMat.block(i,j,h,w).col(c)); for(int r=0; r<h; r++) { - VERIFY_IS_APPROX(m.block(i,j,h,w).col(c).coeff(r), refMat.block(i,j,h,w).col(c).coeff(r)); - } - } - for(int r=0; r<h; r++) - { - VERIFY_IS_APPROX(m.block(i,j,h,w).row(r), refMat.block(i,j,h,w).row(r)); - for(int c=0; c<w; c++) - { - VERIFY_IS_APPROX(m.block(i,j,h,w).row(r).coeff(c), refMat.block(i,j,h,w).row(r).coeff(c)); +// VERIFY_IS_APPROX(m.block(i,j,h,w).col(c).coeff(r), refMat.block(i,j,h,w).col(c).coeff(r)); } } +// for(int r=0; r<h; r++) +// { +// VERIFY_IS_APPROX(m.block(i,j,h,w).row(r), refMat.block(i,j,h,w).row(r)); +// for(int c=0; c<w; c++) +// { +// VERIFY_IS_APPROX(m.block(i,j,h,w).row(r).coeff(c), refMat.block(i,j,h,w).row(r).coeff(c)); +// } +// } } for(int c=0; c<cols; c++) @@ -171,7 +171,7 @@ template<typename Scalar> void sparse_basic(int rows, int cols) } m2.endFill(); std::cerr << m1 << "\n\n" << m2 << "\n"; - VERIFY_IS_APPROX(m1,m2); + VERIFY_IS_APPROX(m2,m1); } // { // m.setZero(); @@ -191,6 +191,17 @@ template<typename Scalar> void sparse_basic(int rows, int cols) // std::cerr << m.transpose() << "\n\n" << refMat.transpose() << "\n\n"; // VERIFY_IS_APPROX(m, refMat); + // test innerVector() + { + DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); + SparseMatrix<Scalar> m2(rows, rows); + initSparse<Scalar>(density, refMat2, m2); + int j0 = ei_random(0,rows-1); + int j1 = ei_random(0,rows-1); +// VERIFY_IS_APPROX(m2.innerVector(j0), refMat2.col(j0)); +// VERIFY_IS_APPROX(m2.innerVector(j0)+m2.innerVector(j1), refMat2.col(j0)+refMat2.col(j1)); + } + // test transpose { DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); diff --git a/test/sparse_vector.cpp b/test/sparse_vector.cpp index 0a25884ca..18ef74833 100644 --- a/test/sparse_vector.cpp +++ b/test/sparse_vector.cpp @@ -62,7 +62,8 @@ template<typename Scalar> void sparse_vector(int rows, int cols) for (typename SparseVectorType::InnerIterator it(v1); it; ++it,++j) { VERIFY(nonzerocoords[j]==it.index()); - VERIFY(it.value()==v1[it.index()]); + VERIFY(it.value()==v1.coeff(it.index())); + VERIFY(it.value()==refV1.coeff(it.index())); } } VERIFY_IS_APPROX(v1, refV1); @@ -76,7 +77,7 @@ template<typename Scalar> void sparse_vector(int rows, int cols) VERIFY_IS_APPROX(v1*s1-v2, refV1*s1-refV2); - std::cerr << v1.dot(v2) << " == " << refV1.dot(refV2) << "\n"; +// std::cerr << v1.dot(v2) << " == " << refV1.dot(refV2) << "\n"; VERIFY_IS_APPROX(v1.dot(v2), refV1.dot(refV2)); } @@ -85,7 +86,7 @@ void test_sparse_vector() { for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST( sparse_vector<double>(8, 8) ); -// CALL_SUBTEST( sparse_vector<std::complex<double> >(16, 16) ); + CALL_SUBTEST( sparse_vector<std::complex<double> >(16, 16) ); CALL_SUBTEST( sparse_vector<double>(299, 535) ); } } |