diff --git a/Eigen/src/SparseCore/SparseBlock.h b/Eigen/src/SparseCore/SparseBlock.h index eefd80702..f13c6dfe2 100644 --- a/Eigen/src/SparseCore/SparseBlock.h +++ b/Eigen/src/SparseCore/SparseBlock.h @@ -12,51 +12,35 @@ namespace Eigen { -namespace internal { -template -struct traits > +template +class BlockImpl + : public SparseMatrixBase > { - typedef typename traits::Scalar Scalar; - typedef typename traits::Index Index; - typedef typename traits::StorageKind StorageKind; - typedef MatrixXpr XprKind; - enum { - IsRowMajor = (int(MatrixType::Flags)&RowMajorBit)==RowMajorBit, - Flags = MatrixType::Flags, - RowsAtCompileTime = IsRowMajor ? Size : MatrixType::RowsAtCompileTime, - ColsAtCompileTime = IsRowMajor ? MatrixType::ColsAtCompileTime : Size, - MaxRowsAtCompileTime = RowsAtCompileTime, - MaxColsAtCompileTime = ColsAtCompileTime, - CoeffReadCost = MatrixType::CoeffReadCost - }; -}; -} // end namespace internal - -template -class SparseInnerVectorSet : internal::no_assignment_operator, - public SparseMatrixBase > -{ - public: - - enum { IsRowMajor = internal::traits::IsRowMajor }; - - EIGEN_SPARSE_PUBLIC_INTERFACE(SparseInnerVectorSet) - class InnerIterator: public MatrixType::InnerIterator + typedef typename internal::remove_all::type _MatrixTypeNested; + typedef Block BlockType; +public: + enum { IsRowMajor = internal::traits::IsRowMajor }; +protected: + enum { OuterSize = IsRowMajor ? BlockRows : BlockCols }; +public: + EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType) + + class InnerIterator: public XprType::InnerIterator { public: - inline InnerIterator(const SparseInnerVectorSet& xpr, Index outer) - : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) + inline InnerIterator(const BlockType& xpr, Index outer) + : XprType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) {} inline Index row() const { return IsRowMajor ? m_outer : this->index(); } inline Index col() const { return IsRowMajor ? this->index() : m_outer; } protected: Index m_outer; }; - class ReverseInnerIterator: public MatrixType::ReverseInnerIterator + class ReverseInnerIterator: public XprType::ReverseInnerIterator { public: - inline ReverseInnerIterator(const SparseInnerVectorSet& xpr, Index outer) - : MatrixType::ReverseInnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) + inline ReverseInnerIterator(const BlockType& xpr, Index outer) + : XprType::ReverseInnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) {} inline Index row() const { return IsRowMajor ? m_outer : this->index(); } inline Index col() const { return IsRowMajor ? this->index() : m_outer; } @@ -64,39 +48,22 @@ class SparseInnerVectorSet : internal::no_assignment_operator, Index m_outer; }; - inline SparseInnerVectorSet(const MatrixType& matrix, Index outerStart, Index outerSize) - : m_matrix(matrix), m_outerStart(outerStart), m_outerSize(outerSize) - { - eigen_assert( (outerStart>=0) && ((outerStart+outerSize)<=matrix.outerSize()) ); - } + inline BlockImpl(const XprType& xpr, int i) + : m_matrix(xpr), m_outerStart(i), m_outerSize(OuterSize) + {} - inline SparseInnerVectorSet(const MatrixType& matrix, Index outer) - : m_matrix(matrix), m_outerStart(outer), m_outerSize(Size) - { - eigen_assert(Size!=Dynamic); - eigen_assert( (outer>=0) && (outer -// inline SparseInnerVectorSet& operator=(const SparseMatrixBase& other) -// { -// return *this; -// } - -// template -// inline SparseInnerVectorSet& operator=(const SparseMatrixBase& other) -// { -// return *this; -// } + inline BlockImpl(const XprType& xpr, int startRow, int startCol, int blockRows, int blockCols) + : m_matrix(xpr), m_outerStart(IsRowMajor ? startRow : startCol), m_outerSize(IsRowMajor ? blockRows : blockCols) + {} EIGEN_STRONG_INLINE Index rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); } EIGEN_STRONG_INLINE Index cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); } protected: - const typename MatrixType::Nested m_matrix; + const typename XprType::Nested m_matrix; Index m_outerStart; - const internal::variable_if_dynamic m_outerSize; + const internal::variable_if_dynamic m_outerSize; }; @@ -104,32 +71,36 @@ class SparseInnerVectorSet : internal::no_assignment_operator, * specialisation for SparseMatrix ***************************************************************************/ -template -class SparseInnerVectorSet, Size> - : public SparseMatrixBase, Size> > +template +class BlockImpl,BlockRows,BlockCols,true,Sparse> + : public SparseMatrixBase,BlockRows,BlockCols,true> > { - typedef SparseMatrix<_Scalar, _Options, _Index> MatrixType; - public: - - enum { IsRowMajor = internal::traits::IsRowMajor }; - - EIGEN_SPARSE_PUBLIC_INTERFACE(SparseInnerVectorSet) - class InnerIterator: public MatrixType::InnerIterator + typedef SparseMatrix<_Scalar, _Options, _Index> XprType; + typedef typename internal::remove_all::type _MatrixTypeNested; + typedef Block BlockType; +public: + enum { IsRowMajor = internal::traits::IsRowMajor }; + EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType) +protected: + enum { OuterSize = IsRowMajor ? BlockRows : BlockCols }; +public: + + class InnerIterator: public XprType::InnerIterator { public: - inline InnerIterator(const SparseInnerVectorSet& xpr, Index outer) - : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) + inline InnerIterator(const BlockType& xpr, Index outer) + : XprType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) {} inline Index row() const { return IsRowMajor ? m_outer : this->index(); } inline Index col() const { return IsRowMajor ? this->index() : m_outer; } protected: Index m_outer; }; - class ReverseInnerIterator: public MatrixType::ReverseInnerIterator + class ReverseInnerIterator: public XprType::ReverseInnerIterator { public: - inline ReverseInnerIterator(const SparseInnerVectorSet& xpr, Index outer) - : MatrixType::ReverseInnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) + inline ReverseInnerIterator(const BlockType& xpr, Index outer) + : XprType::ReverseInnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) {} inline Index row() const { return IsRowMajor ? m_outer : this->index(); } inline Index col() const { return IsRowMajor ? this->index() : m_outer; } @@ -137,23 +108,18 @@ class SparseInnerVectorSet, Size> Index m_outer; }; - inline SparseInnerVectorSet(const MatrixType& matrix, Index outerStart, Index outerSize) - : m_matrix(matrix), m_outerStart(outerStart), m_outerSize(outerSize) - { - eigen_assert( (outerStart>=0) && ((outerStart+outerSize)<=matrix.outerSize()) ); - } + inline BlockImpl(const XprType& xpr, int i) + : m_matrix(xpr), m_outerStart(i), m_outerSize(OuterSize) + {} - inline SparseInnerVectorSet(const MatrixType& matrix, Index outer) - : m_matrix(matrix), m_outerStart(outer), m_outerSize(Size) - { - eigen_assert(Size==1); - eigen_assert( (outer>=0) && (outer - inline SparseInnerVectorSet& operator=(const SparseMatrixBase& other) + inline BlockType& operator=(const SparseMatrixBase& other) { - typedef typename internal::remove_all::type _NestedMatrixType; + typedef typename internal::remove_all::type _NestedMatrixType; _NestedMatrixType& matrix = const_cast<_NestedMatrixType&>(m_matrix);; // This assignement is slow if this vector set is not empty // and/or it is not at the end of the nonzeros of the underlying matrix. @@ -172,7 +138,7 @@ class SparseInnerVectorSet, Size> if(nnz>free_size) { // realloc manually to reduce copies - typename MatrixType::Storage newdata(m_matrix.nonZeros() - nnz_previous + nnz); + typename XprType::Storage newdata(m_matrix.nonZeros() - nnz_previous + nnz); std::memcpy(&newdata.value(0), &m_matrix.data().value(0), nnz_head*sizeof(Scalar)); std::memcpy(&newdata.index(0), &m_matrix.data().index(0), nnz_head*sizeof(Index)); @@ -221,12 +187,12 @@ class SparseInnerVectorSet, Size> matrix.outerIndexPtr()[k] += offset; } - return *this; + return derived(); } - inline SparseInnerVectorSet& operator=(const SparseInnerVectorSet& other) + inline BlockType& operator=(const BlockType& other) { - return operator=(other); + return operator=(other); } inline const Scalar* valuePtr() const @@ -252,12 +218,12 @@ class SparseInnerVectorSet, Size> else if(m_outerSize.value()==0) return 0; else - return Map >(m_matrix.innerNonZeroPtr()+m_outerStart, m_outerSize.value()).sum(); + return Map >(m_matrix.innerNonZeroPtr()+m_outerStart, m_outerSize.value()).sum(); } const Scalar& lastCoeff() const { - EIGEN_STATIC_ASSERT_VECTOR_ONLY(SparseInnerVectorSet); + EIGEN_STATIC_ASSERT_VECTOR_ONLY(BlockImpl); eigen_assert(nonZeros()>0); if(m_matrix.isCompressed()) return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart+1]-1]; @@ -265,122 +231,154 @@ class SparseInnerVectorSet, Size> return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart]+m_matrix.innerNonZeroPtr()[m_outerStart]-1]; } -// template -// inline SparseInnerVectorSet& operator=(const SparseMatrixBase& other) -// { -// return *this; -// } - EIGEN_STRONG_INLINE Index rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); } EIGEN_STRONG_INLINE Index cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); } protected: - typename MatrixType::Nested m_matrix; + typename XprType::Nested m_matrix; Index m_outerStart; - const internal::variable_if_dynamic m_outerSize; + const internal::variable_if_dynamic m_outerSize; }; //---------- -/** \returns the i-th row of the matrix \c *this. For row-major matrix only. */ -template -SparseInnerVectorSet SparseMatrixBase::row(Index i) -{ - EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES); - return innerVector(i); -} - -/** \returns the i-th row of the matrix \c *this. For row-major matrix only. - * (read-only version) */ -template -const SparseInnerVectorSet SparseMatrixBase::row(Index i) const -{ - EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES); - return innerVector(i); -} - -/** \returns the i-th column of the matrix \c *this. For column-major matrix only. */ -template -SparseInnerVectorSet SparseMatrixBase::col(Index i) -{ - EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); - return innerVector(i); -} - -/** \returns the i-th column of the matrix \c *this. For column-major matrix only. - * (read-only version) */ -template -const SparseInnerVectorSet SparseMatrixBase::col(Index i) const -{ - EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); - return innerVector(i); -} - /** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this * is col-major (resp. row-major). */ template -SparseInnerVectorSet SparseMatrixBase::innerVector(Index outer) -{ return SparseInnerVectorSet(derived(), outer); } +typename SparseMatrixBase::InnerVectorReturnType SparseMatrixBase::innerVector(Index outer) +{ return InnerVectorReturnType(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 SparseInnerVectorSet SparseMatrixBase::innerVector(Index outer) const -{ return SparseInnerVectorSet(derived(), outer); } - -/** \returns the i-th row of the matrix \c *this. For row-major matrix only. */ -template -SparseInnerVectorSet SparseMatrixBase::middleRows(Index start, Index size) -{ - EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES); - return innerVectors(start, size); -} - -/** \returns the i-th row of the matrix \c *this. For row-major matrix only. - * (read-only version) */ -template -const SparseInnerVectorSet SparseMatrixBase::middleRows(Index start, Index size) const -{ - EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES); - return innerVectors(start, size); -} - -/** \returns the i-th column of the matrix \c *this. For column-major matrix only. */ -template -SparseInnerVectorSet SparseMatrixBase::middleCols(Index start, Index size) -{ - EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); - return innerVectors(start, size); -} - -/** \returns the i-th column of the matrix \c *this. For column-major matrix only. - * (read-only version) */ -template -const SparseInnerVectorSet SparseMatrixBase::middleCols(Index start, Index size) const -{ - EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); - return innerVectors(start, size); -} - - +const typename SparseMatrixBase::ConstInnerVectorReturnType SparseMatrixBase::innerVector(Index outer) const +{ return ConstInnerVectorReturnType(derived(), outer); } /** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this * is col-major (resp. row-major). */ template -SparseInnerVectorSet SparseMatrixBase::innerVectors(Index outerStart, Index outerSize) -{ return SparseInnerVectorSet(derived(), outerStart, outerSize); } +Block SparseMatrixBase::innerVectors(Index outerStart, Index outerSize) +{ + return Block(derived(), + IsRowMajor ? outerStart : 0, IsRowMajor ? 0 : outerStart, + IsRowMajor ? outerSize : rows(), IsRowMajor ? cols() : outerSize); + +} /** \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 SparseInnerVectorSet SparseMatrixBase::innerVectors(Index outerStart, Index outerSize) const -{ return SparseInnerVectorSet(derived(), outerStart, outerSize); } +const Block SparseMatrixBase::innerVectors(Index outerStart, Index outerSize) const +{ + return Block(derived(), + IsRowMajor ? outerStart : 0, IsRowMajor ? 0 : outerStart, + IsRowMajor ? outerSize : rows(), IsRowMajor ? cols() : outerSize); + +} + +#if 1 + +/** Generic implementation of sparse Block expression. + * Real-only. + */ +template +class BlockImpl + : public SparseMatrixBase >, internal::no_assignment_operator +{ + typedef typename internal::remove_all::type _MatrixTypeNested; + typedef Block BlockType; +public: + enum { IsRowMajor = internal::traits::IsRowMajor }; + EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType) + + /** Column or Row constructor + */ + inline BlockImpl(const XprType& xpr, int i) + : m_matrix(xpr), + m_startRow( (BlockRows==1) && (BlockCols==XprType::ColsAtCompileTime) ? i : 0), + m_startCol( (BlockRows==XprType::RowsAtCompileTime) && (BlockCols==1) ? i : 0), + m_blockRows(xpr.rows()), + m_blockCols(xpr.cols()) + {} + + /** Dynamic-size constructor + */ + inline BlockImpl(const XprType& xpr, int startRow, int startCol, int blockRows, int blockCols) + : m_matrix(xpr), m_startRow(startRow), m_startCol(startCol), m_blockRows(blockRows), m_blockCols(blockCols) + {} + + inline int rows() const { return m_blockRows.value(); } + inline int cols() const { return m_blockCols.value(); } + + inline Scalar& coeffRef(int row, int col) + { + return m_matrix.const_cast_derived() + .coeffRef(row + m_startRow.value(), col + m_startCol.value()); + } + + inline const Scalar coeff(int row, int col) const + { + return m_matrix.coeff(row + m_startRow.value(), col + m_startCol.value()); + } + + inline Scalar& coeffRef(int index) + { + return m_matrix.const_cast_derived() + .coeffRef(m_startRow.value() + (RowsAtCompileTime == 1 ? 0 : index), + m_startCol.value() + (RowsAtCompileTime == 1 ? index : 0)); + } + + inline const Scalar coeff(int index) const + { + return m_matrix + .coeff(m_startRow.value() + (RowsAtCompileTime == 1 ? 0 : index), + m_startCol.value() + (RowsAtCompileTime == 1 ? index : 0)); + } + + inline const _MatrixTypeNested& nestedExpression() const { return m_matrix; } + + class InnerIterator : public _MatrixTypeNested::InnerIterator + { + typedef typename _MatrixTypeNested::InnerIterator Base; + const BlockType& m_block; + Index m_end; + public: + + EIGEN_STRONG_INLINE InnerIterator(const BlockType& block, Index outer) + : Base(block.derived().nestedExpression(), outer + (IsRowMajor ? block.m_startRow.value() : block.m_startCol.value())), + m_block(block), + m_end(IsRowMajor ? block.m_startCol.value()+block.m_blockCols.value() : block.m_startRow.value()+block.m_blockRows.value()) + { + while(Base::index() < (IsRowMajor ? m_block.m_startCol.value() : m_block.m_startRow.value())) + Base::operator++(); + } + + inline Index index() const { return Base::index() - (IsRowMajor ? m_block.m_startCol.value() : m_block.m_startRow.value()); } + inline Index outer() const { return Base::outer() - (IsRowMajor ? m_block.m_startRow.value() : m_block.m_startCol.value()); } + inline Index row() const { return Base::row() - m_block.m_startRow.value(); } + inline Index col() const { return Base::col() - m_block.m_startCol.value(); } + + inline operator bool() const { return Base::operator bool() && Base::index() < m_end; } + }; + + protected: + friend class InnerIterator; + + const typename XprType::Nested m_matrix; + const internal::variable_if_dynamic m_startRow; + const internal::variable_if_dynamic m_startCol; + const internal::variable_if_dynamic m_blockRows; + const internal::variable_if_dynamic m_blockCols; + +}; + +#endif } // end namespace Eigen diff --git a/Eigen/src/SparseCore/SparseDiagonalProduct.h b/Eigen/src/SparseCore/SparseDiagonalProduct.h index ccba02124..5ec4018e6 100644 --- a/Eigen/src/SparseCore/SparseDiagonalProduct.h +++ b/Eigen/src/SparseCore/SparseDiagonalProduct.h @@ -118,12 +118,12 @@ class sparse_diagonal_product_inner_iterator_selector : public CwiseBinaryOp< scalar_product_op, - SparseInnerVectorSet, + typename Rhs::ConstInnerVectorReturnType, typename Lhs::DiagonalVectorType>::InnerIterator { typedef typename CwiseBinaryOp< scalar_product_op, - SparseInnerVectorSet, + typename Rhs::ConstInnerVectorReturnType, typename Lhs::DiagonalVectorType>::InnerIterator Base; typedef typename Lhs::Index Index; Index m_outer; @@ -156,12 +156,12 @@ class sparse_diagonal_product_inner_iterator_selector : public CwiseBinaryOp< scalar_product_op, - SparseInnerVectorSet, + typename Lhs::ConstInnerVectorReturnType, Transpose >::InnerIterator { typedef typename CwiseBinaryOp< scalar_product_op, - SparseInnerVectorSet, + typename Lhs::ConstInnerVectorReturnType, Transpose >::InnerIterator Base; typedef typename Lhs::Index Index; Index m_outer; diff --git a/Eigen/src/SparseCore/SparseMatrixBase.h b/Eigen/src/SparseCore/SparseMatrixBase.h index 39332a178..2ef841616 100644 --- a/Eigen/src/SparseCore/SparseMatrixBase.h +++ b/Eigen/src/SparseCore/SparseMatrixBase.h @@ -136,13 +136,13 @@ template class SparseMatrixBase : public EigenBase # include "../plugins/CommonCwiseBinaryOps.h" # include "../plugins/MatrixCwiseUnaryOps.h" # include "../plugins/MatrixCwiseBinaryOps.h" +# include "../plugins/BlockMethods.h" # ifdef EIGEN_SPARSEMATRIXBASE_PLUGIN # include EIGEN_SPARSEMATRIXBASE_PLUGIN # endif # undef EIGEN_CURRENT_STORAGE_BASE_CLASS #undef EIGEN_CURRENT_STORAGE_BASE_CLASS - /** \returns the number of rows. \sa cols() */ inline Index rows() const { return derived().rows(); } /** \returns the number of columns. \sa rows() */ @@ -392,26 +392,15 @@ template class SparseMatrixBase : public EigenBase const Transpose transpose() const { return derived(); } const AdjointReturnType adjoint() const { return transpose(); } - // sub-vector - SparseInnerVectorSet row(Index i); - const SparseInnerVectorSet row(Index i) const; - SparseInnerVectorSet col(Index j); - const SparseInnerVectorSet col(Index j) const; - SparseInnerVectorSet innerVector(Index outer); - const SparseInnerVectorSet innerVector(Index outer) const; + // inner-vector + typedef Block InnerVectorReturnType; + typedef Block ConstInnerVectorReturnType; + InnerVectorReturnType innerVector(Index outer); + const ConstInnerVectorReturnType innerVector(Index outer) const; - // set of sub-vectors - SparseInnerVectorSet subrows(Index start, Index size); - const SparseInnerVectorSet subrows(Index start, Index size) const; - SparseInnerVectorSet subcols(Index start, Index size); - const SparseInnerVectorSet subcols(Index start, Index size) const; - - SparseInnerVectorSet middleRows(Index start, Index size); - const SparseInnerVectorSet middleRows(Index start, Index size) const; - SparseInnerVectorSet middleCols(Index start, Index size); - const SparseInnerVectorSet middleCols(Index start, Index size) const; - SparseInnerVectorSet innerVectors(Index outerStart, Index outerSize); - const SparseInnerVectorSet innerVectors(Index outerStart, Index outerSize) const; + // set of inner-vectors + Block innerVectors(Index outerStart, Index outerSize); + const Block innerVectors(Index outerStart, Index outerSize) const; /** \internal use operator= */ template diff --git a/Eigen/src/SparseCore/SparseUtil.h b/Eigen/src/SparseCore/SparseUtil.h index a686e08da..d58b51356 100644 --- a/Eigen/src/SparseCore/SparseUtil.h +++ b/Eigen/src/SparseCore/SparseUtil.h @@ -73,7 +73,6 @@ template class Dynamic template class SparseVector; template class MappedSparseMatrix; -template class SparseInnerVectorSet; template class SparseTriangularView; template class SparseSelfAdjointView; template class SparseDiagonalProduct; diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp index 4566de9f2..cebb5e6a0 100644 --- a/test/sparse_basic.cpp +++ b/test/sparse_basic.cpp @@ -155,63 +155,6 @@ template void sparse_basic(const SparseMatrixType& re VERIFY_IS_APPROX(m2,m1); } - // test basic computations - { - DenseMatrix refM1 = DenseMatrix::Zero(rows, rows); - DenseMatrix refM2 = DenseMatrix::Zero(rows, rows); - DenseMatrix refM3 = DenseMatrix::Zero(rows, rows); - DenseMatrix refM4 = DenseMatrix::Zero(rows, rows); - SparseMatrixType m1(rows, rows); - SparseMatrixType m2(rows, rows); - SparseMatrixType m3(rows, rows); - SparseMatrixType m4(rows, rows); - initSparse(density, refM1, m1); - initSparse(density, refM2, m2); - initSparse(density, refM3, m3); - initSparse(density, refM4, m4); - - VERIFY_IS_APPROX(m1+m2, refM1+refM2); - VERIFY_IS_APPROX(m1+m2+m3, refM1+refM2+refM3); - VERIFY_IS_APPROX(m3.cwiseProduct(m1+m2), refM3.cwiseProduct(refM1+refM2)); - VERIFY_IS_APPROX(m1*s1-m2, refM1*s1-refM2); - - VERIFY_IS_APPROX(m1*=s1, refM1*=s1); - VERIFY_IS_APPROX(m1/=s1, refM1/=s1); - - VERIFY_IS_APPROX(m1+=m2, refM1+=refM2); - VERIFY_IS_APPROX(m1-=m2, refM1-=refM2); - - if(SparseMatrixType::IsRowMajor) - VERIFY_IS_APPROX(m1.innerVector(0).dot(refM2.row(0)), refM1.row(0).dot(refM2.row(0))); - else - VERIFY_IS_APPROX(m1.innerVector(0).dot(refM2.row(0)), refM1.col(0).dot(refM2.row(0))); - - VERIFY_IS_APPROX(m1.conjugate(), refM1.conjugate()); - VERIFY_IS_APPROX(m1.real(), refM1.real()); - - refM4.setRandom(); - // sparse cwise* dense - VERIFY_IS_APPROX(m3.cwiseProduct(refM4), refM3.cwiseProduct(refM4)); -// VERIFY_IS_APPROX(m3.cwise()/refM4, refM3.cwise()/refM4); - - // test aliasing - VERIFY_IS_APPROX((m1 = -m1), (refM1 = -refM1)); - VERIFY_IS_APPROX((m1 = m1.transpose()), (refM1 = refM1.transpose().eval())); - VERIFY_IS_APPROX((m1 = -m1.transpose()), (refM1 = -refM1.transpose().eval())); - VERIFY_IS_APPROX((m1 += -m1), (refM1 += -refM1)); - } - - // test transpose - { - DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); - SparseMatrixType m2(rows, rows); - initSparse(density, refMat2, m2); - VERIFY_IS_APPROX(m2.transpose().eval(), refMat2.transpose().eval()); - VERIFY_IS_APPROX(m2.transpose(), refMat2.transpose()); - - VERIFY_IS_APPROX(SparseMatrixType(m2.adjoint()), refMat2.adjoint()); - } - // test innerVector() { DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); @@ -274,6 +217,86 @@ template void sparse_basic(const SparseMatrixType& re //m2.innerVectors(j0,n0) = m2.innerVectors(j0,n0) + m2.innerVectors(j1,n0); //refMat2.block(0,j0,rows,n0) = refMat2.block(0,j0,rows,n0) + refMat2.block(0,j1,rows,n0); } + + // test basic computations + { + DenseMatrix refM1 = DenseMatrix::Zero(rows, rows); + DenseMatrix refM2 = DenseMatrix::Zero(rows, rows); + DenseMatrix refM3 = DenseMatrix::Zero(rows, rows); + DenseMatrix refM4 = DenseMatrix::Zero(rows, rows); + SparseMatrixType m1(rows, rows); + SparseMatrixType m2(rows, rows); + SparseMatrixType m3(rows, rows); + SparseMatrixType m4(rows, rows); + initSparse(density, refM1, m1); + initSparse(density, refM2, m2); + initSparse(density, refM3, m3); + initSparse(density, refM4, m4); + + VERIFY_IS_APPROX(m1+m2, refM1+refM2); + VERIFY_IS_APPROX(m1+m2+m3, refM1+refM2+refM3); + VERIFY_IS_APPROX(m3.cwiseProduct(m1+m2), refM3.cwiseProduct(refM1+refM2)); + VERIFY_IS_APPROX(m1*s1-m2, refM1*s1-refM2); + + VERIFY_IS_APPROX(m1*=s1, refM1*=s1); + VERIFY_IS_APPROX(m1/=s1, refM1/=s1); + + VERIFY_IS_APPROX(m1+=m2, refM1+=refM2); + VERIFY_IS_APPROX(m1-=m2, refM1-=refM2); + + if(SparseMatrixType::IsRowMajor) + VERIFY_IS_APPROX(m1.innerVector(0).dot(refM2.row(0)), refM1.row(0).dot(refM2.row(0))); + else + VERIFY_IS_APPROX(m1.innerVector(0).dot(refM2.row(0)), refM1.col(0).dot(refM2.row(0))); + + VERIFY_IS_APPROX(m1.conjugate(), refM1.conjugate()); + VERIFY_IS_APPROX(m1.real(), refM1.real()); + + refM4.setRandom(); + // sparse cwise* dense + VERIFY_IS_APPROX(m3.cwiseProduct(refM4), refM3.cwiseProduct(refM4)); +// VERIFY_IS_APPROX(m3.cwise()/refM4, refM3.cwise()/refM4); + + // test aliasing + VERIFY_IS_APPROX((m1 = -m1), (refM1 = -refM1)); + VERIFY_IS_APPROX((m1 = m1.transpose()), (refM1 = refM1.transpose().eval())); + VERIFY_IS_APPROX((m1 = -m1.transpose()), (refM1 = -refM1.transpose().eval())); + VERIFY_IS_APPROX((m1 += -m1), (refM1 += -refM1)); + } + + // test transpose + { + DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); + SparseMatrixType m2(rows, rows); + initSparse(density, refMat2, m2); + VERIFY_IS_APPROX(m2.transpose().eval(), refMat2.transpose().eval()); + VERIFY_IS_APPROX(m2.transpose(), refMat2.transpose()); + + VERIFY_IS_APPROX(SparseMatrixType(m2.adjoint()), refMat2.adjoint()); + } + + + + // test generic blocks + { + DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); + SparseMatrixType m2(rows, rows); + initSparse(density, refMat2, m2); + int j0 = internal::random(0,rows-2); + int j1 = internal::random(0,rows-2); + int n0 = internal::random(1,rows-(std::max)(j0,j1)); + if(SparseMatrixType::IsRowMajor) + VERIFY_IS_APPROX(m2.block(j0,0,n0,cols), refMat2.block(j0,0,n0,cols)); + else + VERIFY_IS_APPROX(m2.block(0,j0,rows,n0), refMat2.block(0,j0,rows,n0)); + + if(SparseMatrixType::IsRowMajor) + VERIFY_IS_APPROX(m2.block(j0,0,n0,cols)+m2.block(j1,0,n0,cols), + refMat2.block(j0,0,n0,cols)+refMat2.block(j1,0,n0,cols)); + else + VERIFY_IS_APPROX(m2.block(0,j0,rows,n0)+m2.block(0,j1,rows,n0), + refMat2.block(0,j0,rows,n0)+refMat2.block(0,j1,rows,n0)); + } // test prune {