From c4c70669d165afefe0c68e7bb194ee81b9fba0b5 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 14 Jan 2009 14:24:10 +0000 Subject: [PATCH] Big rewrite in the Sparse module: SparseMatrixBase no longer inherits MatrixBase. That means a lot of features which were available for sparse matrices via the dense (and super slow) implemention are no longer available. All features which make sense for sparse matrices (aka can be implemented efficiently) will be implemented soon, but don't expect to see an API as rich as for the dense path. Other changes: * no block(), row(), col() anymore. * instead use .innerVector() to get a col or row vector of a matrix. * .segment(), start(), end() will be back soon, not sure for block() * faster cwise product --- CMakeLists.txt | 6 + Eigen/Sparse | 8 + Eigen/src/Core/CwiseBinaryOp.h | 4 +- Eigen/src/Core/CwiseUnaryOp.h | 2 - Eigen/src/Core/Dot.h | 11 +- Eigen/src/Core/MatrixBase.h | 4 +- Eigen/src/Core/Product.h | 3 - Eigen/src/Core/Sum.h | 11 +- Eigen/src/Core/Transpose.h | 4 +- Eigen/src/Core/util/Constants.h | 2 +- Eigen/src/Sparse/CoreIterators.h | 82 +--- Eigen/src/Sparse/SparseAssign.h | 0 Eigen/src/Sparse/SparseBlock.h | 75 +++- Eigen/src/Sparse/SparseCwise.h | 181 ++++++++ Eigen/src/Sparse/SparseCwiseBinaryOp.h | 560 +++++++++++++++++++++++++ Eigen/src/Sparse/SparseCwiseUnaryOp.h | 175 ++++++++ Eigen/src/Sparse/SparseDot.h | 97 +++++ Eigen/src/Sparse/SparseFlagged.h | 98 +++++ Eigen/src/Sparse/SparseFuzzy.h | 41 ++ Eigen/src/Sparse/SparseMatrix.h | 8 +- Eigen/src/Sparse/SparseMatrixBase.h | 412 +++++++++++++++++- Eigen/src/Sparse/SparseProduct.h | 70 ++-- Eigen/src/Sparse/SparseRedux.h | 50 +-- Eigen/src/Sparse/SparseTranspose.h | 85 ++++ Eigen/src/Sparse/SparseUtil.h | 50 +++ Eigen/src/Sparse/SparseVector.h | 2 +- Eigen/src/Sparse/TriangularSolver.h | 31 ++ cmake/FindSuperLU.cmake | 4 - test/CMakeLists.txt | 2 +- test/sparse_basic.cpp | 35 +- test/sparse_vector.cpp | 7 +- 31 files changed, 1920 insertions(+), 200 deletions(-) create mode 100644 Eigen/src/Sparse/SparseAssign.h create mode 100644 Eigen/src/Sparse/SparseCwise.h create mode 100644 Eigen/src/Sparse/SparseCwiseBinaryOp.h create mode 100644 Eigen/src/Sparse/SparseCwiseUnaryOp.h create mode 100644 Eigen/src/Sparse/SparseDot.h create mode 100644 Eigen/src/Sparse/SparseFlagged.h create mode 100644 Eigen/src/Sparse/SparseFuzzy.h create mode 100644 Eigen/src/Sparse/SparseTranspose.h 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::LhsNested LhsNested; typedef typename ei_traits::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 template::Vectorization, - int Unrolling = ei_dot_traits::Unrolling, - int Storage = (ei_traits::Flags | ei_traits::Flags) & SparseBit + int Unrolling = ei_dot_traits::Unrolling > struct ei_dot_impl; template -struct ei_dot_impl +struct ei_dot_impl { typedef typename Derived1::Scalar Scalar; static Scalar run(const Derived1& v1, const Derived2& v2) @@ -164,12 +163,12 @@ struct ei_dot_impl }; template -struct ei_dot_impl +struct ei_dot_impl : public ei_dot_novec_unroller {}; template -struct ei_dot_impl +struct ei_dot_impl { typedef typename Derived1::Scalar Scalar; typedef typename ei_packet_traits::type PacketScalar; @@ -222,7 +221,7 @@ struct ei_dot_impl -struct ei_dot_impl +struct ei_dot_impl { typedef typename Derived1::Scalar Scalar; typedef typename ei_packet_traits::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 class MatrixBase { return lazyAssign(other._expression()); } /** Overloaded for sparse product evaluation */ - template - Derived& lazyAssign(const Product& product); + /*template + Derived& lazyAssign(const Product& product);*/ CommaInitializer 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 * - NormalProduct * - CacheFriendlyProduct * - DiagonalProduct - * - SparseProduct */ template struct ei_product_mode { @@ -87,8 +86,6 @@ template 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 template::Vectorization, - int Unrolling = ei_sum_traits::Unrolling, - int Storage = ei_traits::Flags & SparseBit + int Unrolling = ei_sum_traits::Unrolling > struct ei_sum_impl; template -struct ei_sum_impl +struct ei_sum_impl { typedef typename Derived::Scalar Scalar; static Scalar run(const Derived& mat) @@ -178,12 +177,12 @@ struct ei_sum_impl }; template -struct ei_sum_impl +struct ei_sum_impl : public ei_sum_novec_unroller {}; template -struct ei_sum_impl +struct ei_sum_impl { typedef typename Derived::Scalar Scalar; typedef typename ei_packet_traits::type PacketScalar; @@ -228,7 +227,7 @@ struct ei_sum_impl }; template -struct ei_sum_impl +struct ei_sum_impl { typedef typename Derived::Scalar Scalar; typedef typename ei_packet_traits::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 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 { // 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 -class MatrixBase::InnerIterator +/** \class InnerIterator + * \brief An InnerIterator allows to loop over the element of a sparse (or dense) matrix or expression + * + * todo + */ +// template InnerIterator; + +// generic version for dense matrix and expressions +template class MatrixBase::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::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 -class Transpose::InnerIterator : public MatrixType::InnerIterator -{ - public: - - EIGEN_STRONG_INLINE InnerIterator(const Transpose& trans, int outer) - : MatrixType::InnerIterator(trans.m_matrix, outer) - {} -}; +/* template class Block::InnerIterator { @@ -138,50 +137,9 @@ class Block::InnerIter int m_start; int m_end; int m_offset; -}; - -template -class CwiseUnaryOp::InnerIterator -{ - typedef typename CwiseUnaryOp::Scalar Scalar; - typedef typename ei_traits::_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 struct ei_is_scalar_product { enum { ret = false }; }; template struct ei_is_scalar_product > { enum { ret = true }; }; @@ -204,8 +162,8 @@ class CwiseBinaryOp::InnerIterator EIGEN_STRONG_INLINE InnerIterator(const CwiseBinaryOp& binOp, int outer) : Base(binOp.m_lhs,binOp.m_rhs,binOp.m_functor,outer) {} -}; - +};*/ +/* template class CwiseBinaryOpInnerIterator { @@ -270,7 +228,7 @@ class CwiseBinaryOpInnerIterator Scalar m_value; int m_id; }; -/* + template class CwiseBinaryOpInnerIterator,Lhs,Rhs,Derived> { diff --git a/Eigen/src/Sparse/SparseAssign.h b/Eigen/src/Sparse/SparseAssign.h new file mode 100644 index 000000000..e69de29bb 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 . -#ifndef EIGEN_SPARSEBLOCK_H -#define EIGEN_SPARSEBLOCK_H +#ifndef EIGEN_SPARSE_BLOCK_H +#define EIGEN_SPARSE_BLOCK_H +template +struct ei_traits > +{ + typedef typename ei_traits::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 +class SparseInnerVector : ei_no_assignment_operator, + public SparseMatrixBase > +{ + enum { + IsRowMajor = ei_traits::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 +class SparseInnerVector::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 +SparseInnerVector SparseMatrixBase::innerVector(int outer) +{ return SparseInnerVector(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 +const SparseInnerVector SparseMatrixBase::innerVector(int outer) const +{ return SparseInnerVector(derived(), outer); } + +# if 0 template class Block : public SparseMatrixBase > @@ -117,6 +184,6 @@ public: const ei_int_if_dynamic 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 +// Copyright (C) 2008 Benoit Jacob +// +// 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 . + +#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::Scalar>, ExpressionType, OtherDerived> + +#define EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE \ + SparseCwiseBinaryOp< \ + ei_scalar_product_op< \ + typename ei_scalar_product_traits< \ + typename ei_traits::Scalar, \ + typename ei_traits::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::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::Scalar>, ExpressionType, \ + NestByValue >*/ + +template class SparseCwise +{ + public: + + typedef typename ei_traits::Scalar Scalar; + typedef typename ei_meta_if::ret, + ExpressionType, const ExpressionType&>::ret ExpressionTypeNested; + typedef CwiseUnaryOp, ExpressionType> ScalarAddReturnType; + + inline SparseCwise(const ExpressionType& matrix) : m_matrix(matrix) {} + + /** \internal */ + inline const ExpressionType& _expression() const { return m_matrix; } + + template + const EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE + operator*(const SparseMatrixBase &other) const; + + template + const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op) + operator/(const SparseMatrixBase &other) const; + + template + const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_min_op) + min(const SparseMatrixBase &other) const; + + template + const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_max_op) + max(const SparseMatrixBase &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 + inline ExpressionType& operator*=(const MatrixBase &other); + + template + inline ExpressionType& operator/=(const MatrixBase &other); + + template const EIGEN_CWISE_BINOP_RETURN_TYPE(std::less) + operator<(const MatrixBase& other) const; + + template const EIGEN_CWISE_BINOP_RETURN_TYPE(std::less_equal) + operator<=(const MatrixBase& other) const; + + template const EIGEN_CWISE_BINOP_RETURN_TYPE(std::greater) + operator>(const MatrixBase& other) const; + + template const EIGEN_CWISE_BINOP_RETURN_TYPE(std::greater_equal) + operator>=(const MatrixBase& other) const; + + template const EIGEN_CWISE_BINOP_RETURN_TYPE(std::equal_to) + operator==(const MatrixBase& other) const; + + template const EIGEN_CWISE_BINOP_RETURN_TYPE(std::not_equal_to) + operator!=(const MatrixBase& 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 +inline const SparseCwise +SparseMatrixBase::cwise() const +{ + return derived(); +} + +template +inline SparseCwise +SparseMatrixBase::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 +// +// 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 . + +#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 class ei_sparse_cwise_binary_op +// : ei_no_assignment_operator +// { +// typedef CwiseBinaryOp CwiseBinaryXpr; +// public: +// +// typedef typename ei_traits::LhsNested LhsNested; +// typedef typename ei_traits::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::ret +// ? int(ei_is_same_type::ret) +// : int(ei_is_same_type::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 +struct ei_traits > +{ + 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::type _LhsNested; + typedef typename ei_unref::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::Cost + }; +}; + +template +class SparseCwiseBinaryOp : ei_no_assignment_operator, + public SparseMatrixBase > +{ + public: + + class InnerIterator; + + EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseCwiseBinaryOp) + typedef typename ei_traits::LhsNested LhsNested; + typedef typename ei_traits::RhsNested RhsNested; + typedef typename ei_unref::type _LhsNested; + typedef typename ei_unref::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::ret + ? int(ei_is_same_type::ret) + : int(ei_is_same_type::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 +class CwiseBinaryOp, Rhs> + : public ei_sparse_cwise_binary_op, + public SparseMatrixBase, 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(lhs, rhs, func) + {} +}; +template +class CwiseBinaryOp, SparseMatrixBase > + : public ei_sparse_cwise_binary_op, + public SparseMatrixBase > +{ + 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(lhs, rhs, func) + {} +};*/ + +// template +// class CwiseBinaryOp > : ei_no_assignment_operator, +// public SparseMatrixBase, SparseMatrixBase > > +// { +// public: +// +// EIGEN_GENERIC_PUBLIC_INTERFACE(CwiseBinaryOp) +// typedef typename ei_traits::LhsNested LhsNested; +// typedef typename ei_traits::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::ret +// ? int(ei_is_same_type::ret) +// : int(ei_is_same_type::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 +// class CwiseBinaryOp, SparseMatrixBase > : ei_no_assignment_operator, +// public SparseMatrixBase, SparseMatrixBase > > +// { +// public: +// +// EIGEN_GENERIC_PUBLIC_INTERFACE(CwiseBinaryOp) +// typedef typename ei_traits::LhsNested LhsNested; +// typedef typename ei_traits::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::ret +// ? int(ei_is_same_type::ret) +// : int(ei_is_same_type::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 struct ei_is_scalar_product { enum { ret = false }; }; +// template struct ei_is_scalar_product > { enum { ret = true }; }; + +template +class ei_sparse_cwise_binary_op_inner_iterator_selector; + +template +class SparseCwiseBinaryOp::InnerIterator + : public ei_sparse_cwise_binary_op_inner_iterator_selector::InnerIterator> +{ + typedef ei_sparse_cwise_binary_op_inner_iterator_selector< + BinaryOp,Lhs,Rhs, typename SparseCwiseBinaryOp::InnerIterator> Base; + public: + typedef typename SparseCwiseBinaryOp::Scalar Scalar; + typedef typename ei_traits::_LhsNested _LhsNested; + typedef typename _LhsNested::InnerIterator LhsIterator; + typedef typename ei_traits::_RhsNested _RhsNested; + typedef typename _RhsNested::InnerIterator RhsIterator; + + EIGEN_STRONG_INLINE InnerIterator(const SparseCwiseBinaryOp& binOp, int outer) + : Base(binOp,outer) + {} +}; + +/*************************************************************************** +* Implementation of inner-iterators +***************************************************************************/ + +// template struct ei_is_scalar_product { enum { ret = false }; }; +// template struct ei_is_scalar_product > { enum { ret = true }; }; + +// helper class + + +// sparse - sparse (generic) +template +class ei_sparse_cwise_binary_op_inner_iterator_selector +{ + typedef SparseCwiseBinaryOp CwiseBinaryXpr; + typedef typename ei_traits::Scalar Scalar; + typedef typename ei_traits::_LhsNested _LhsNested; + typedef typename ei_traits::_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(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 +// class SparseCwiseBinaryOp::InnerIterator +// : public ei_sparse_cwise_binary_op_inner_iterator_selector::InnerIterator> +// { +// typedef CwiseBinaryOpInnerIterator< +// BinaryOp,Lhs,Rhs, typename CwiseBinaryOp::InnerIterator> Base; +// public: +// typedef typename CwiseBinaryOp::Scalar Scalar; +// typedef typename ei_traits::_LhsNested _LhsNested; +// typedef typename _LhsNested::InnerIterator LhsIterator; +// typedef typename ei_traits::_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 +class ei_sparse_cwise_binary_op_inner_iterator_selector, Lhs, Rhs, Derived, IsSparse, IsSparse> +{ + typedef ei_scalar_product_op BinaryOp; + typedef SparseCwiseBinaryOp CwiseBinaryXpr; + typedef typename CwiseBinaryXpr::Scalar Scalar; + typedef typename ei_traits::_LhsNested _LhsNested; + typedef typename _LhsNested::InnerIterator LhsIterator; + typedef typename ei_traits::_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(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 +class ei_sparse_cwise_binary_op_inner_iterator_selector, Lhs, Rhs, Derived, IsSparse, IsDense> +{ + typedef ei_scalar_product_op BinaryOp; + typedef SparseCwiseBinaryOp CwiseBinaryXpr; + typedef typename CwiseBinaryXpr::Scalar Scalar; + typedef typename ei_traits::_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(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 +class ei_sparse_cwise_binary_op_inner_iterator_selector, Lhs, Rhs, Derived, IsDense, IsSparse> +{ + typedef ei_scalar_product_op BinaryOp; + typedef SparseCwiseBinaryOp CwiseBinaryXpr; + typedef typename CwiseBinaryXpr::Scalar Scalar; + typedef typename ei_traits::_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(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 +template +EIGEN_STRONG_INLINE const SparseCwiseBinaryOp::Scalar>, + Derived, OtherDerived> +SparseMatrixBase::operator-(const SparseMatrixBase &other) const +{ + return SparseCwiseBinaryOp, + Derived, OtherDerived>(derived(), other.derived()); +} + +template +template +EIGEN_STRONG_INLINE Derived & +SparseMatrixBase::operator-=(const SparseMatrixBase &other) +{ + return *this = *this - other; +} + +template +template +EIGEN_STRONG_INLINE const SparseCwiseBinaryOp::Scalar>, Derived, OtherDerived> +SparseMatrixBase::operator+(const SparseMatrixBase &other) const +{ + return SparseCwiseBinaryOp, Derived, OtherDerived>(derived(), other.derived()); +} + +template +template +EIGEN_STRONG_INLINE Derived & +SparseMatrixBase::operator+=(const SparseMatrixBase& other) +{ + return *this = *this + other; +} + +template +template +EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE +SparseCwise::operator*(const SparseMatrixBase &other) const +{ + return EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE(_expression(), other.derived()); +} + +template +template +EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op) +SparseCwise::operator/(const SparseMatrixBase &other) const +{ + return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op)(_expression(), other.derived()); +} + +// template +// template +// inline ExpressionType& Cwise::operator*=(const SparseMatrixBase &other) +// { +// return m_matrix.const_cast_derived() = *this * other; +// } + +// template +// template +// inline ExpressionType& Cwise::operator/=(const SparseMatrixBase &other) +// { +// return m_matrix.const_cast_derived() = *this / other; +// } + +template +template +EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_min_op) +SparseCwise::min(const SparseMatrixBase &other) const +{ + return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_min_op)(_expression(), other.derived()); +} + +template +template +EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_max_op) +SparseCwise::max(const SparseMatrixBase &other) const +{ + return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_max_op)(_expression(), other.derived()); +} + +// template +// template +// EIGEN_STRONG_INLINE const CwiseBinaryOp +// SparseMatrixBase::binaryExpr(const SparseMatrixBase &other, const CustomBinaryOp& func) const +// { +// return CwiseBinaryOp(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 +// +// 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 . + +#ifndef EIGEN_SPARSE_CWISE_UNARY_OP_H +#define EIGEN_SPARSE_CWISE_UNARY_OP_H + +template +struct ei_traits > : ei_traits +{ + typedef typename ei_result_of< + UnaryOp(typename MatrixType::Scalar) + >::type Scalar; + typedef typename MatrixType::Nested MatrixTypeNested; + typedef typename ei_unref::type _MatrixTypeNested; + enum { + CoeffReadCost = _MatrixTypeNested::CoeffReadCost + ei_functor_traits::Cost + }; +}; + +template +class SparseCwiseUnaryOp : ei_no_assignment_operator, + public SparseMatrixBase > +{ + public: + + class InnerIterator; +// typedef typename ei_unref::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 +class SparseCwiseUnaryOp::InnerIterator +{ + typedef typename SparseCwiseUnaryOp::Scalar Scalar; + typedef typename ei_traits::_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 +template +EIGEN_STRONG_INLINE const SparseCwiseUnaryOp +SparseMatrixBase::unaryExpr(const CustomUnaryOp& func) const +{ + return SparseCwiseUnaryOp(derived(), func); +} + +template +EIGEN_STRONG_INLINE const SparseCwiseUnaryOp::Scalar>,Derived> +SparseMatrixBase::operator-() const +{ + return derived(); +} + +template +EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_abs_op) +SparseCwise::abs() const +{ + return _expression(); +} + +template +EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_abs2_op) +SparseCwise::abs2() const +{ + return _expression(); +} + +template +EIGEN_STRONG_INLINE typename SparseMatrixBase::ConjugateReturnType +SparseMatrixBase::conjugate() const +{ + return ConjugateReturnType(derived()); +} + +template +EIGEN_STRONG_INLINE const typename SparseMatrixBase::RealReturnType +SparseMatrixBase::real() const { return derived(); } + +template +EIGEN_STRONG_INLINE const typename SparseMatrixBase::ImagReturnType +SparseMatrixBase::imag() const { return derived(); } + +template +template +EIGEN_STRONG_INLINE const SparseCwiseUnaryOp::Scalar, NewType>, Derived> +SparseMatrixBase::cast() const +{ + return derived(); +} + +template +EIGEN_STRONG_INLINE const SparseCwiseUnaryOp::Scalar>, Derived> +SparseMatrixBase::operator*(const Scalar& scalar) const +{ + return SparseCwiseUnaryOp, Derived> + (derived(), ei_scalar_multiple_op(scalar)); +} + +template +EIGEN_STRONG_INLINE const SparseCwiseUnaryOp::Scalar>, Derived> +SparseMatrixBase::operator/(const Scalar& scalar) const +{ + return SparseCwiseUnaryOp, Derived> + (derived(), ei_scalar_quotient1_op(scalar)); +} + +template +EIGEN_STRONG_INLINE Derived& +SparseMatrixBase::operator*=(const Scalar& other) +{ + return *this = *this * other; +} + +template +EIGEN_STRONG_INLINE Derived& +SparseMatrixBase::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 +// +// 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 . + +#ifndef EIGEN_SPARSE_DOT_H +#define EIGEN_SPARSE_DOT_H + +template +template +typename ei_traits::Scalar +SparseMatrixBase::dot(const MatrixBase& 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::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 +template +typename ei_traits::Scalar +SparseMatrixBase::dot(const SparseMatrixBase& 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::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() +inline typename NumTraits::Scalar>::Real +SparseMatrixBase::squaredNorm() const +{ + return ei_real((*this).cwise().abs2().sum()); +} + +template +inline typename NumTraits::Scalar>::Real +SparseMatrixBase::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 +// Copyright (C) 2008 Gael Guennebaud +// +// 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 . + +#ifndef EIGEN_SPARSE_FLAGGED_H +#define EIGEN_SPARSE_FLAGGED_H + +template +struct ei_traits > : ei_traits +{ + enum { Flags = (ExpressionType::Flags | Added) & ~Removed }; +}; + +template class SparseFlagged + : public SparseMatrixBase > +{ + public: + + EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseFlagged) + class InnerIterator; + class ReverseInnerIterator; + + typedef typename ei_meta_if::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 + class SparseFlagged::InnerIterator : public ExpressionType::InnerIterator +{ + public: + + EIGEN_STRONG_INLINE InnerIterator(const SparseFlagged& xpr, int outer) + : ExpressionType::InnerIterator(xpr.m_matrix, outer) + {} +}; + +template + class SparseFlagged::ReverseInnerIterator : public ExpressionType::ReverseInnerIterator +{ + public: + + EIGEN_STRONG_INLINE ReverseInnerIterator(const SparseFlagged& xpr, int outer) + : ExpressionType::ReverseInnerIterator(xpr.m_matrix, outer) + {} +}; + +template +template +inline const SparseFlagged +SparseMatrixBase::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 +// +// 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 . + +#ifndef EIGEN_SPARSE_FUZZY_H +#define EIGEN_SPARSE_FUZZY_H + +// template +// template +// bool SparseMatrixBase::isApprox( +// const OtherDerived& other, +// typename NumTraits::Real prec +// ) const +// { +// const typename ei_nested::type nested(derived()); +// const typename ei_nested::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 > { public: - EIGEN_GENERIC_PUBLIC_INTERFACE(SparseMatrix) + EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseMatrix) protected: public: typedef SparseMatrixBase SparseBase; enum { - RowMajor = SparseBase::RowMajor + RowMajor = SparseBase::IsRowMajor }; typedef SparseMatrix TransposedSparseMatrix; @@ -267,7 +267,7 @@ class SparseMatrix } template - inline SparseMatrix(const MatrixBase& other) + inline SparseMatrix(const SparseMatrixBase& other) : m_outerSize(0), m_innerSize(0), m_outerIndex(0) { *this = other.derived(); @@ -305,7 +305,7 @@ class SparseMatrix } template - inline SparseMatrix& operator=(const MatrixBase& other) + inline SparseMatrix& operator=(const SparseMatrixBase& 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 -class SparseMatrixBase : public MatrixBase +template class SparseMatrixBase { public: - typedef MatrixBase Base; - typedef typename Base::Scalar Scalar; - typedef typename Base::RealScalar RealScalar; + typedef typename ei_traits::Scalar Scalar; +// typedef typename Derived::InnerIterator InnerIterator; + enum { - Flags = Base::Flags, - RowMajor = ei_traits::Flags&RowMajorBit ? 1 : 0 + + RowsAtCompileTime = ei_traits::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::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::RowsAtCompileTime, + ei_traits::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::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::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::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::IsComplex, + const SparseCwiseUnaryOp, Derived>, + const Derived& + >::ret ConjugateReturnType; + /** \internal the return type of MatrixBase::real() */ + typedef CwiseUnaryOp, Derived> RealReturnType; + /** \internal the return type of MatrixBase::imag() */ + typedef CwiseUnaryOp, Derived> ImagReturnType; + /** \internal the return type of MatrixBase::adjoint() */ + typedef Eigen::Transpose::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 then RealScalar is \a T. + * + * \sa class NumTraits + */ + typedef typename NumTraits::Real RealScalar; + + /** type of the equivalent square matrix */ + typedef Matrix SquareMatrixType; + inline const Derived& derived() const { return *static_cast(this); } inline Derived& derived() { return *static_cast(this); } inline Derived& const_cast_derived() const { return *static_cast(const_cast(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=(other); return derived(); } + template - inline Derived& operator=(const MatrixBase& other) + inline void assignGeneric(const OtherDerived& other) { // std::cout << "Derived& operator=(const MatrixBase& other)\n"; //const bool transpose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit); @@ -87,9 +175,9 @@ class SparseMatrixBase : public MatrixBase temp.endFill(); derived() = temp.markAsRValue(); - return derived(); } + template inline Derived& operator=(const SparseMatrixBase& other) { @@ -110,7 +198,7 @@ class SparseMatrixBase : public MatrixBase 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().endFill(); } else - this->operator=(static_cast&>(other)); + { + assignGeneric(other.derived()); + } return derived(); } template - inline Derived& operator=(const Product& product); + inline Derived& operator=(const SparseProduct& product); friend std::ostream & operator << (std::ostream & s, const SparseMatrixBase& m) { @@ -167,6 +257,292 @@ class SparseMatrixBase : public MatrixBase return s; } + const SparseCwiseUnaryOp::Scalar>,Derived> operator-() const; + + template + const SparseCwiseBinaryOp::Scalar>, Derived, OtherDerived> + operator+(const SparseMatrixBase &other) const; + + template + const SparseCwiseBinaryOp::Scalar>, Derived, OtherDerived> + operator-(const SparseMatrixBase &other) const; + + template + Derived& operator+=(const SparseMatrixBase& other); + template + Derived& operator-=(const SparseMatrixBase& other); + +// template +// Derived& operator+=(const Flagged, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit>& other); + + Derived& operator*=(const Scalar& other); + Derived& operator/=(const Scalar& other); + + const SparseCwiseUnaryOp::Scalar>, Derived> + operator*(const Scalar& scalar) const; + const SparseCwiseUnaryOp::Scalar>, Derived> + operator/(const Scalar& scalar) const; + + inline friend const SparseCwiseUnaryOp::Scalar>, Derived> + operator*(const Scalar& scalar, const SparseMatrixBase& matrix) + { return matrix*scalar; } + + + template + const typename SparseProductReturnType::Type + operator*(const SparseMatrixBase &other) const; + + template + Derived& operator*=(const SparseMatrixBase& other); + + template + typename ei_plain_matrix_type_column_major::type + solveTriangular(const MatrixBase& other) const; + + template + void solveTriangularInPlace(MatrixBase& other) const; + + template Scalar dot(const MatrixBase& other) const; + template Scalar dot(const SparseMatrixBase& other) const; + RealScalar squaredNorm() const; + RealScalar norm() const; +// const PlainMatrixType normalized() const; +// void normalize(); + + SparseTranspose transpose() { return derived(); } + const SparseTranspose transpose() const { return derived(); } + // void transposeInPlace(); + // const AdjointReturnType adjoint() const; + + SparseInnerVector innerVector(int outer); + const SparseInnerVector 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::Type block(int startRow, int startCol, int blockRows, int blockCols); +// const typename BlockReturnType::Type +// block(int startRow, int startCol, int blockRows, int blockCols) const; +// +// typename BlockReturnType::SubVectorType segment(int start, int size); +// const typename BlockReturnType::SubVectorType segment(int start, int size) const; +// +// typename BlockReturnType::SubVectorType start(int size); +// const typename BlockReturnType::SubVectorType start(int size) const; +// +// typename BlockReturnType::SubVectorType end(int size); +// const typename BlockReturnType::SubVectorType end(int size) const; +// +// typename BlockReturnType::Type corner(CornerType type, int cRows, int cCols); +// const typename BlockReturnType::Type corner(CornerType type, int cRows, int cCols) const; +// +// template +// typename BlockReturnType::Type block(int startRow, int startCol); +// template +// const typename BlockReturnType::Type block(int startRow, int startCol) const; + +// template +// typename BlockReturnType::Type corner(CornerType type); +// template +// const typename BlockReturnType::Type corner(CornerType type) const; + +// template typename BlockReturnType::SubVectorType start(void); +// template const typename BlockReturnType::SubVectorType start() const; + +// template typename BlockReturnType::SubVectorType end(); +// template const typename BlockReturnType::SubVectorType end() const; + +// template typename BlockReturnType::SubVectorType segment(int start); +// template const typename BlockReturnType::SubVectorType segment(int start) const; + +// DiagonalCoeffs diagonal(); +// const DiagonalCoeffs diagonal() const; + +// template Part part(); +// template const Part 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 +// static const CwiseNullaryOp NullaryExpr(int rows, int cols, const CustomNullaryOp& func); +// template +// static const CwiseNullaryOp NullaryExpr(int size, const CustomNullaryOp& func); +// template +// static const CwiseNullaryOp 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 asDiagonal() const; + +// Derived& setConstant(const Scalar& value); +// Derived& setZero(); +// Derived& setOnes(); +// Derived& setRandom(); +// Derived& setIdentity(); + + Matrix toDense() const + { + Matrix res(rows(),cols()); + res.setZero(); + for (int j=0; j + bool isApprox(const SparseMatrixBase& other, + RealScalar prec = precision()) const + { return toDense().isApprox(other.toDense(),prec); } + + template + bool isApprox(const MatrixBase& other, + RealScalar prec = precision()) const + { return toDense().isApprox(other,prec); } +// bool isMuchSmallerThan(const RealScalar& other, +// RealScalar prec = precision()) const; +// template +// bool isMuchSmallerThan(const MatrixBase& other, +// RealScalar prec = precision()) const; + +// bool isApproxToConstant(const Scalar& value, RealScalar prec = precision()) const; +// bool isZero(RealScalar prec = precision()) const; +// bool isOnes(RealScalar prec = precision()) const; +// bool isIdentity(RealScalar prec = precision()) const; +// bool isDiagonal(RealScalar prec = precision()) const; + +// bool isUpperTriangular(RealScalar prec = precision()) const; +// bool isLowerTriangular(RealScalar prec = precision()) const; + +// template +// bool isOrthogonal(const MatrixBase& other, +// RealScalar prec = precision()) const; +// bool isUnitary(RealScalar prec = precision()) const; + +// template +// inline bool operator==(const MatrixBase& other) const +// { return (cwise() == other).all(); } + +// template +// inline bool operator!=(const MatrixBase& other) const +// { return (cwise() != other).any(); } + + + template + const SparseCwiseUnaryOp::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::type eval() const + { return typename ei_eval::type(derived()); } + +// template +// void swap(const MatrixBase& other); + + template + const SparseFlagged marked() const; +// const Flagged 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 nestByValue() const; + + + ConjugateReturnType conjugate() const; + const RealReturnType real() const; + const ImagReturnType imag() const; + + template + const SparseCwiseUnaryOp unaryExpr(const CustomUnaryOp& func = CustomUnaryOp()) const; + +// template +// const CwiseBinaryOp +// binaryExpr(const MatrixBase &other, const CustomBinaryOp& func = CustomBinaryOp()) const; + + + Scalar sum() const; +// Scalar trace() const; + +// typename ei_traits::Scalar minCoeff() const; +// typename ei_traits::Scalar maxCoeff() const; + +// typename ei_traits::Scalar minCoeff(int* row, int* col = 0) const; +// typename ei_traits::Scalar maxCoeff(int* row, int* col = 0) const; + +// template +// typename ei_result_of::Scalar)>::type +// redux(const BinaryOp& func) const; + +// template +// void visit(Visitor& func) const; + + + const SparseCwise cwise() const; + SparseCwise cwise(); + +// inline const WithFormat format(const IOFormat& fmt) const; + +/////////// Array module /////////// + /* + bool all(void) const; + bool any(void) const; + + const PartialRedux rowwise() const; + const PartialRedux colwise() const; + + static const CwiseNullaryOp,Derived> Random(int rows, int cols); + static const CwiseNullaryOp,Derived> Random(int size); + static const CwiseNullaryOp,Derived> Random(); + + template + const Select + select(const MatrixBase& thenMatrix, + const MatrixBase& elseMatrix) const; + + template + inline const Select > + select(const MatrixBase& thenMatrix, typename ThenDerived::Scalar elseScalar) const; + + template + inline const Select, ElseDerived > + select(typename ElseDerived::Scalar thenScalar, const MatrixBase& elseMatrix) const; + + template RealScalar lpNorm() const; + */ + + // template // Scalar dot(const MatrixBase& 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 -struct ProductReturnType +struct SparseProductReturnType { typedef typename ei_traits::Scalar Scalar; enum { @@ -47,12 +47,11 @@ struct ProductReturnType SparseMatrix, const typename ei_nested::type>::ret RhsNested; - typedef Product Type; + typedef SparseProduct Type; }; template -struct ei_traits > +struct ei_traits > { // clean the nested types: typedef typename ei_cleantype::type _LhsNested; @@ -87,30 +86,28 @@ struct ei_traits > }; }; -template class Product : ei_no_assignment_operator, - public MatrixBase > +template +class SparseProduct : ei_no_assignment_operator, + public SparseMatrixBase > { public: - EIGEN_GENERIC_PUBLIC_INTERFACE(Product) + EIGEN_GENERIC_PUBLIC_INTERFACE(SparseProduct) private: - typedef typename ei_traits::_LhsNested _LhsNested; - typedef typename ei_traits::_RhsNested _RhsNested; + typedef typename ei_traits::_LhsNested _LhsNested; + typedef typename ei_traits::_RhsNested _RhsNested; public: template - 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,21 +228,21 @@ struct ei_sparse_product_selector // by ProductReturnType which transform it to (col col *) by evaluating rhs. -template -template -inline Derived& MatrixBase::lazyAssign(const Product& product) -{ -// std::cout << "sparse product to dense\n"; - ei_sparse_product_selector< - typename ei_cleantype::type, - typename ei_cleantype::type, - typename ei_cleantype::type>::run(product.lhs(),product.rhs(),derived()); - return derived(); -} +// template +// template +// inline Derived& SparseMatrixBase::lazyAssign(const SparseProduct& product) +// { +// // std::cout << "sparse product to dense\n"; +// ei_sparse_product_selector< +// typename ei_cleantype::type, +// typename ei_cleantype::type, +// typename ei_cleantype::type>::run(product.lhs(),product.rhs(),derived()); +// return derived(); +// } template template -inline Derived& SparseMatrixBase::operator=(const Product& product) +inline Derived& SparseMatrixBase::operator=(const SparseProduct& product) { // std::cout << "sparse product to sparse\n"; ei_sparse_product_selector< @@ -255,4 +252,27 @@ inline Derived& SparseMatrixBase::operator=(const Product +template +inline const typename SparseProductReturnType::Type +SparseMatrixBase::operator*(const SparseMatrixBase &other) const +{ + 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::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 -struct ei_sum_impl +template +typename ei_traits::Scalar +SparseMatrixBase::sum() const { - 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 -struct ei_dot_impl -{ - 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()0 && cols()>0 && "you are using a non initialized matrix"); + Scalar res = 0; + for (int j=0; j +// +// 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 . + +#ifndef EIGEN_SPARSETRANSPOSE_H +#define EIGEN_SPARSETRANSPOSE_H + +template +struct ei_traits > : ei_traits > +{}; + +template class SparseTranspose + : public SparseMatrixBase > +{ + 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 class SparseTranspose::InnerIterator : public MatrixType::InnerIterator +{ + public: + + EIGEN_STRONG_INLINE InnerIterator(const SparseTranspose& trans, int outer) + : MatrixType::InnerIterator(trans.m_matrix, outer) + {} +}; + +template class SparseTranspose::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 \ +EIGEN_STRONG_INLINE Derived& operator Op(const Eigen::SparseMatrixBase& 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 \ +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::Scalar Scalar; \ +typedef typename Eigen::NumTraits::Real RealScalar; \ +typedef typename Eigen::ei_nested::type Nested; \ +enum { RowsAtCompileTime = Eigen::ei_traits::RowsAtCompileTime, \ + ColsAtCompileTime = Eigen::ei_traits::ColsAtCompileTime, \ + Flags = Eigen::ei_traits::Flags, \ + CoeffReadCost = Eigen::ei_traits::CoeffReadCost, \ + SizeAtCompileTime = Base::SizeAtCompileTime, \ + IsVectorAtCompileTime = Base::IsVectorAtCompileTime }; + +#define EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(Derived) \ +_EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(Derived, Eigen::SparseMatrixBase) + enum SparseBackend { DefaultBackend, Taucs, @@ -64,6 +104,16 @@ template class SparseMatrixBase; template class SparseMatrix; template class SparseVector; +template class SparseTranspose; +template class SparseInnerVector; +template class SparseCwise; +template class SparseCwiseUnaryOp; +template class SparseCwiseBinaryOp; +template class SparseProduct; +template class SparseFlagged; + +template 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 > { 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 } }; +template +template +void SparseMatrixBase::solveTriangularInPlace(MatrixBase& 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::Flags & RowMajorBit }; + + typedef typename ei_meta_if::type, OtherDerived&>::ret OtherCopy; + OtherCopy otherCopy(other.derived()); + + ei_solve_triangular_selector::type>::run(derived(), otherCopy); + + if (copy) + other = otherCopy; +} + +template +template +typename ei_plain_matrix_type_column_major::type +SparseMatrixBase::solveTriangular(const MatrixBase& other) const +{ + typename ei_plain_matrix_type_column_major::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 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 void sparse_basic(int rows, int cols) int w = ei_random(1,cols-j-1); int h = ei_random(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 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 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 m2(rows, rows); + initSparse(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 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 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(8, 8) ); -// CALL_SUBTEST( sparse_vector >(16, 16) ); + CALL_SUBTEST( sparse_vector >(16, 16) ); CALL_SUBTEST( sparse_vector(299, 535) ); } }