From a9688f0b717568713ff1acef3466a0da00a68980 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 9 Feb 2009 09:59:30 +0000 Subject: [PATCH] - add diagonal * sparse product as an expression - split sparse_basic unit test - various fixes in sparse module --- Eigen/Sparse | 1 + Eigen/src/Sparse/CholmodSupport.h | 2 +- Eigen/src/Sparse/DynamicSparseMatrix.h | 2 +- Eigen/src/Sparse/SparseBlock.h | 97 +++++++++++++- Eigen/src/Sparse/SparseCwiseBinaryOp.h | 20 +-- Eigen/src/Sparse/SparseCwiseUnaryOp.h | 2 +- Eigen/src/Sparse/SparseDiagonalProduct.h | 157 +++++++++++++++++++++++ Eigen/src/Sparse/SparseMatrix.h | 2 +- Eigen/src/Sparse/SparseProduct.h | 25 +++- Eigen/src/Sparse/SparseUtil.h | 1 + Eigen/src/Sparse/SparseVector.h | 2 +- Eigen/src/Sparse/SuperLUSupport.h | 2 +- Eigen/src/Sparse/TaucsSupport.h | 2 +- Eigen/src/Sparse/UmfPackSupport.h | 2 +- test/CMakeLists.txt | 1 + test/sparse_basic.cpp | 65 +--------- test/sparse_product.cpp | 132 +++++++++++++++++++ test/sparse_vector.cpp | 1 + 18 files changed, 423 insertions(+), 93 deletions(-) create mode 100644 Eigen/src/Sparse/SparseDiagonalProduct.h create mode 100644 test/sparse_product.cpp diff --git a/Eigen/Sparse b/Eigen/Sparse index 289985ed9..536c28454 100644 --- a/Eigen/Sparse +++ b/Eigen/Sparse @@ -103,6 +103,7 @@ namespace Eigen { #include "src/Sparse/SparseFuzzy.h" #include "src/Sparse/SparseFlagged.h" #include "src/Sparse/SparseProduct.h" +#include "src/Sparse/SparseDiagonalProduct.h" #include "src/Sparse/TriangularSolver.h" #include "src/Sparse/SparseLLT.h" #include "src/Sparse/SparseLDLT.h" diff --git a/Eigen/src/Sparse/CholmodSupport.h b/Eigen/src/Sparse/CholmodSupport.h index 1550fd055..dfd9c787a 100644 --- a/Eigen/src/Sparse/CholmodSupport.h +++ b/Eigen/src/Sparse/CholmodSupport.h @@ -1,7 +1,7 @@ // 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-2009 Gael Guennebaud // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public diff --git a/Eigen/src/Sparse/DynamicSparseMatrix.h b/Eigen/src/Sparse/DynamicSparseMatrix.h index 74916e0ab..72930bfcc 100644 --- a/Eigen/src/Sparse/DynamicSparseMatrix.h +++ b/Eigen/src/Sparse/DynamicSparseMatrix.h @@ -1,7 +1,7 @@ // 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-2009 Gael Guennebaud // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public diff --git a/Eigen/src/Sparse/SparseBlock.h b/Eigen/src/Sparse/SparseBlock.h index a90ccd812..e50852553 100644 --- a/Eigen/src/Sparse/SparseBlock.h +++ b/Eigen/src/Sparse/SparseBlock.h @@ -91,9 +91,9 @@ class SparseInnerVectorSet : ei_no_assignment_operator, }; -//---------- -// specialisation for DynamicSparseMatrix -//---------- +/*************************************************************************** +* specialisation for DynamicSparseMatrix +***************************************************************************/ template class SparseInnerVectorSet, Size> @@ -169,6 +169,89 @@ class SparseInnerVectorSet, Size> }; +/*************************************************************************** +* specialisation for SparseMatrix +***************************************************************************/ +/* +template +class SparseInnerVectorSet, Size> + : public SparseMatrixBase, Size> > +{ + typedef DynamicSparseMatrix<_Scalar, _Options> MatrixType; + enum { IsRowMajor = ei_traits::IsRowMajor }; + public: + + EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseInnerVectorSet) + class InnerIterator: public MatrixType::InnerIterator + { + public: + inline InnerIterator(const SparseInnerVectorSet& xpr, int outer) + : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer) + {} + }; + + inline SparseInnerVectorSet(const MatrixType& matrix, int outerStart, int outerSize) + : m_matrix(matrix), m_outerStart(outerStart), m_outerSize(outerSize) + { + ei_assert( (outerStart>=0) && ((outerStart+outerSize)<=matrix.outerSize()) ); + } + + inline SparseInnerVectorSet(const MatrixType& matrix, int outer) + : m_matrix(matrix), m_outerStart(outer) + { + ei_assert(Size==1); + ei_assert( (outer>=0) && (outer + inline SparseInnerVectorSet& operator=(const SparseMatrixBase& other) + { + if (IsRowMajor != ((OtherDerived::Flags&RowMajorBit)==RowMajorBit)) + { + // need to transpose => perform a block evaluation followed by a big swap + DynamicSparseMatrix aux(other); + *this = aux.markAsRValue(); + } + else + { + // evaluate/copy vector per vector + for (int j=0; j aux(other.innerVector(j)); + m_matrix.const_cast_derived()._data()[m_outerStart+j].swap(aux._data()); + } + } + return *this; + } + + inline SparseInnerVectorSet& operator=(const SparseInnerVectorSet& other) + { + return operator=(other); + } + + inline const Scalar* _valuePtr() const + { return m_matrix._valuePtr() + m_matrix._outerIndexPtr()[m_outerStart]; } + inline const int* _innerIndexPtr() const + { return m_matrix._innerIndexPtr() + m_matrix._outerIndexPtr()[m_outerStart]; } + inline const int* _outerIndexPtr() const { return m_matrix._outerIndexPtr() + m_outerStart; } + +// template +// inline SparseInnerVectorSet& operator=(const SparseMatrixBase& other) +// { +// return *this; +// } + + EIGEN_STRONG_INLINE int rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); } + EIGEN_STRONG_INLINE int cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); } + + protected: + + const typename MatrixType::Nested m_matrix; + int m_outerStart; + const ei_int_if_dynamic m_outerSize; + +}; +*/ //---------- /** \returns the i-th row of the matrix \c *this. For row-major matrix only. */ @@ -192,7 +275,7 @@ const SparseInnerVectorSet SparseMatrixBase::row(int i) cons template SparseInnerVectorSet SparseMatrixBase::col(int i) { - EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COL_MAJOR_MATRICES); + EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); return innerVector(i); } @@ -201,7 +284,7 @@ SparseInnerVectorSet SparseMatrixBase::col(int i) template const SparseInnerVectorSet SparseMatrixBase::col(int i) const { - EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COL_MAJOR_MATRICES); + EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); return innerVector(i); } @@ -242,7 +325,7 @@ const SparseInnerVectorSet SparseMatrixBase::subrows(i template SparseInnerVectorSet SparseMatrixBase::subcols(int start, int size) { - EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COL_MAJOR_MATRICES); + EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); return innerVectors(start, size); } @@ -251,7 +334,7 @@ SparseInnerVectorSet SparseMatrixBase::subcols(int sta template const SparseInnerVectorSet SparseMatrixBase::subcols(int start, int size) const { - EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COL_MAJOR_MATRICES); + EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); return innerVectors(start, size); } diff --git a/Eigen/src/Sparse/SparseCwiseBinaryOp.h b/Eigen/src/Sparse/SparseCwiseBinaryOp.h index a061368ab..d19970efc 100644 --- a/Eigen/src/Sparse/SparseCwiseBinaryOp.h +++ b/Eigen/src/Sparse/SparseCwiseBinaryOp.h @@ -132,11 +132,10 @@ class SparseCwiseBinaryOp::InnerIterator * Implementation of inner-iterators ***************************************************************************/ -// template struct ei_is_scalar_product { enum { ret = false }; }; -// template struct ei_is_scalar_product > { enum { ret = true }; }; - -// helper class +// template struct ei_func_is_conjunction { enum { ret = false }; }; +// template struct ei_func_is_conjunction > { enum { ret = true }; }; +// TODO generalize the ei_scalar_product_op specialization to all conjunctions if any ! // sparse - sparse (generic) template @@ -261,12 +260,13 @@ class ei_sparse_cwise_binary_op_inner_iterator_selector, typedef SparseCwiseBinaryOp CwiseBinaryXpr; typedef typename CwiseBinaryXpr::Scalar Scalar; typedef typename ei_traits::_LhsNested _LhsNested; + typedef typename ei_traits::RhsNested RhsNested; 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) + : m_rhs(xpr.rhs()), m_lhsIter(xpr.lhs(),outer), m_functor(xpr.functor()), m_outer(outer) {} EIGEN_STRONG_INLINE Derived& operator++() @@ -277,7 +277,7 @@ class ei_sparse_cwise_binary_op_inner_iterator_selector, 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)); } + m_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 int row() const { return m_lhsIter.row(); } @@ -286,9 +286,9 @@ class ei_sparse_cwise_binary_op_inner_iterator_selector, EIGEN_STRONG_INLINE operator bool() const { return m_lhsIter; } protected: - const CwiseBinaryXpr& m_xpr; + const RhsNested m_rhs; LhsIterator m_lhsIter; - const BinaryFunc& m_functor; + const BinaryFunc m_functor; const int m_outer; }; @@ -331,6 +331,10 @@ class ei_sparse_cwise_binary_op_inner_iterator_selector, }; +/*************************************************************************** +* Implementation of SparseMatrixBase and SparseCwise functions/operators +***************************************************************************/ + template template EIGEN_STRONG_INLINE const SparseCwiseBinaryOp::Scalar>, diff --git a/Eigen/src/Sparse/SparseCwiseUnaryOp.h b/Eigen/src/Sparse/SparseCwiseUnaryOp.h index ff7c9c73f..b11c0f8a3 100644 --- a/Eigen/src/Sparse/SparseCwiseUnaryOp.h +++ b/Eigen/src/Sparse/SparseCwiseUnaryOp.h @@ -89,7 +89,7 @@ class SparseCwiseUnaryOp::InnerIterator protected: MatrixTypeIterator m_iter; - const UnaryOp& m_functor; + const UnaryOp m_functor; }; template diff --git a/Eigen/src/Sparse/SparseDiagonalProduct.h b/Eigen/src/Sparse/SparseDiagonalProduct.h new file mode 100644 index 000000000..983823a0c --- /dev/null +++ b/Eigen/src/Sparse/SparseDiagonalProduct.h @@ -0,0 +1,157 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2009 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_DIAGONAL_PRODUCT_H +#define EIGEN_SPARSE_DIAGONAL_PRODUCT_H + +// the product a diagonal matrix with a sparse matrix can be easily +// implemented using expression template. We have two very different cases: +// 1 - diag * row-major sparse +// => each inner vector <=> scalar * sparse vector product +// => so we can reuse CwiseUnaryOp::InnerIterator +// 2 - diag * col-major sparse +// => each inner vector <=> densevector * sparse vector cwise product +// => again, we can reuse specialization of CwiseBinaryOp::InnerIterator +// for that particular case +// The two other cases are symmetric. + +template +struct ei_traits > : ei_traits > +{ + typedef typename ei_cleantype::type _Lhs; + typedef typename ei_cleantype::type _Rhs; + enum { + SparseFlags = ((int(_Lhs::Flags)&Diagonal)==Diagonal) ? _Rhs::Flags : _Lhs::Flags, + Flags = SparseBit | (SparseFlags&RowMajorBit) + }; +}; + +enum {SDP_IsDiagonal, SDP_IsSparseRowMajor, SDP_IsSparseColMajor}; +template +class ei_sparse_diagonal_product_inner_iterator_selector; + +template +class SparseDiagonalProduct : public SparseMatrixBase >, ei_no_assignment_operator +{ + typedef typename ei_traits::_LhsNested _LhsNested; + typedef typename ei_traits::_RhsNested _RhsNested; + + enum { + LhsMode = (_LhsNested::Flags&Diagonal)==Diagonal ? SDP_IsDiagonal + : (_LhsNested::Flags&RowMajorBit) ? SDP_IsSparseRowMajor : SDP_IsSparseColMajor, + RhsMode = (_RhsNested::Flags&Diagonal)==Diagonal ? SDP_IsDiagonal + : (_RhsNested::Flags&RowMajorBit) ? SDP_IsSparseRowMajor : SDP_IsSparseColMajor + }; + + public: + + EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseDiagonalProduct) + + typedef ei_sparse_diagonal_product_inner_iterator_selector + <_LhsNested,_RhsNested,SparseDiagonalProduct,LhsMode,RhsMode> InnerIterator; + + template + EIGEN_STRONG_INLINE SparseDiagonalProduct(const Lhs& lhs, const Rhs& rhs) + : m_lhs(lhs), m_rhs(rhs) + { + ei_assert(lhs.cols() == rhs.rows() && "invalid sparse matrix * diagonal matrix product"); + } + + EIGEN_STRONG_INLINE int rows() const { return m_lhs.rows(); } + EIGEN_STRONG_INLINE int cols() const { return m_rhs.cols(); } + + EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; } + EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; } + + protected: + LhsNested m_lhs; + RhsNested m_rhs; +}; + + +template +class ei_sparse_diagonal_product_inner_iterator_selector + + : public SparseCwiseUnaryOp,Rhs>::InnerIterator +{ + typedef typename SparseCwiseUnaryOp,Rhs>::InnerIterator Base; + public: + inline ei_sparse_diagonal_product_inner_iterator_selector( + const SparseDiagonalProductType& expr, int outer) + : Base(expr.rhs()*(expr.lhs().diagonal().coeff(outer)), outer) + {} +}; + +template +class ei_sparse_diagonal_product_inner_iterator_selector + + : public SparseCwiseBinaryOp< + ei_scalar_product_op, + SparseInnerVectorSet, + typename Lhs::_CoeffsVectorType>::InnerIterator +{ + typedef typename SparseCwiseBinaryOp< + ei_scalar_product_op, + SparseInnerVectorSet, + typename Lhs::_CoeffsVectorType>::InnerIterator Base; + public: + inline ei_sparse_diagonal_product_inner_iterator_selector( + const SparseDiagonalProductType& expr, int outer) + : Base(expr.rhs().innerVector(outer) .cwise()* expr.lhs().diagonal(), 0) + {} +}; + +template +class ei_sparse_diagonal_product_inner_iterator_selector + + : public SparseCwiseUnaryOp,Lhs>::InnerIterator +{ + typedef typename SparseCwiseUnaryOp,Lhs>::InnerIterator Base; + public: + inline ei_sparse_diagonal_product_inner_iterator_selector( + const SparseDiagonalProductType& expr, int outer) + : Base(expr.lhs()*expr.rhs().diagonal().coeff(outer), outer) + {} +}; + +template +class ei_sparse_diagonal_product_inner_iterator_selector + + : public SparseCwiseBinaryOp< + ei_scalar_product_op, + SparseInnerVectorSet, + NestByValue > >::InnerIterator +{ + typedef typename SparseCwiseBinaryOp< + ei_scalar_product_op, + SparseInnerVectorSet, + NestByValue > >::InnerIterator Base; + public: + inline ei_sparse_diagonal_product_inner_iterator_selector( + const SparseDiagonalProductType& expr, int outer) + : Base(expr.lhs().innerVector(outer) .cwise()* expr.rhs().diagonal().transpose().nestByValue(), 0) + {} +}; + +#endif // EIGEN_SPARSE_DIAGONAL_PRODUCT_H diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h index 23766d516..8aebe4ae2 100644 --- a/Eigen/src/Sparse/SparseMatrix.h +++ b/Eigen/src/Sparse/SparseMatrix.h @@ -1,7 +1,7 @@ // 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-2009 Gael Guennebaud // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public diff --git a/Eigen/src/Sparse/SparseProduct.h b/Eigen/src/Sparse/SparseProduct.h index f3cf465f1..d3ef57455 100644 --- a/Eigen/src/Sparse/SparseProduct.h +++ b/Eigen/src/Sparse/SparseProduct.h @@ -29,7 +29,9 @@ template struct ei_sparse_product_mode { enum { - value = (Rhs::Flags&Lhs::Flags&SparseBit)==SparseBit + value = ((Lhs::Flags&Diagonal)==Diagonal || (Rhs::Flags&Diagonal)==Diagonal) + ? DiagonalProduct + : (Rhs::Flags&Lhs::Flags&SparseBit)==SparseBit ? SparseTimeSparseProduct : (Lhs::Flags&SparseBit)==SparseBit ? SparseTimeDenseProduct @@ -45,6 +47,15 @@ struct SparseProductReturnType typedef SparseProduct Type; }; +template +struct SparseProductReturnType +{ + typedef const typename ei_nested::type LhsNested; + typedef const typename ei_nested::type RhsNested; + + typedef SparseDiagonalProduct Type; +}; + // sparse product return type specialization template struct SparseProductReturnType @@ -95,7 +106,7 @@ struct ei_traits > // RhsIsRowMajor = (RhsFlags & RowMajorBit)==RowMajorBit, EvalToRowMajor = (RhsFlags & LhsFlags & RowMajorBit), - ResultIsSparse = ProductMode==SparseTimeSparseProduct, + ResultIsSparse = ProductMode==SparseTimeSparseProduct || ProductMode==DiagonalProduct, RemovedBits = ~( (EvalToRowMajor ? 0 : RowMajorBit) | (ResultIsSparse ? 0 : SparseBit) ), @@ -112,7 +123,8 @@ struct ei_traits > }; template -class SparseProduct : ei_no_assignment_operator, public ei_traits >::Base +class SparseProduct : ei_no_assignment_operator, + public ei_traits >::Base { public: @@ -285,7 +297,6 @@ template template inline Derived& SparseMatrixBase::operator=(const SparseProduct& product) { -// std::cout << "sparse product to sparse\n"; ei_sparse_product_selector< typename ei_cleantype::type, typename ei_cleantype::type, @@ -311,7 +322,7 @@ inline Derived& SparseMatrixBase::operator=(const SparseProduct template Derived& MatrixBase::lazyAssign(const SparseProduct& product) -{ +{ typedef typename ei_cleantype::type _Lhs; typedef typename ei_cleantype::type _Rhs; typedef typename _Lhs::InnerIterator LhsInnerIterator; @@ -333,7 +344,7 @@ Derived& MatrixBase::lazyAssign(const SparseProduct foo = derived().row(j); + Block res(derived().row(LhsIsRowMajor ? j : 0)); for (; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i) { if (LhsIsSelfAdjoint) @@ -345,7 +356,7 @@ Derived& MatrixBase::lazyAssign(const SparseProduct class SparseCwiseUnaryO template class SparseCwiseBinaryOp; template class SparseFlagged; +template class SparseDiagonalProduct; template struct ei_sparse_product_mode; template::value> struct SparseProductReturnType; diff --git a/Eigen/src/Sparse/SparseVector.h b/Eigen/src/Sparse/SparseVector.h index a5100150b..8e5a6efed 100644 --- a/Eigen/src/Sparse/SparseVector.h +++ b/Eigen/src/Sparse/SparseVector.h @@ -1,7 +1,7 @@ // 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-2009 Gael Guennebaud // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public diff --git a/Eigen/src/Sparse/SuperLUSupport.h b/Eigen/src/Sparse/SuperLUSupport.h index ded145eeb..825de0bb5 100644 --- a/Eigen/src/Sparse/SuperLUSupport.h +++ b/Eigen/src/Sparse/SuperLUSupport.h @@ -1,7 +1,7 @@ // 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-2009 Gael Guennebaud // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public diff --git a/Eigen/src/Sparse/TaucsSupport.h b/Eigen/src/Sparse/TaucsSupport.h index 2035a2af7..cda26a529 100644 --- a/Eigen/src/Sparse/TaucsSupport.h +++ b/Eigen/src/Sparse/TaucsSupport.h @@ -1,7 +1,7 @@ // 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-2009 Gael Guennebaud // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public diff --git a/Eigen/src/Sparse/UmfPackSupport.h b/Eigen/src/Sparse/UmfPackSupport.h index 1bb40a7c3..b76ffb252 100644 --- a/Eigen/src/Sparse/UmfPackSupport.h +++ b/Eigen/src/Sparse/UmfPackSupport.h @@ -1,7 +1,7 @@ // 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-2009 Gael Guennebaud // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 99a21835f..4d2ac7501 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -146,6 +146,7 @@ if(QT4_FOUND) endif(QT4_FOUND) ei_add_test(sparse_vector) ei_add_test(sparse_basic) +ei_add_test(sparse_product) ei_add_test(sparse_solvers " " "${SPARSE_LIBS}") ei_add_test(reverse) diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp index 23b8526b7..410ef96a6 100644 --- a/test/sparse_basic.cpp +++ b/test/sparse_basic.cpp @@ -238,6 +238,8 @@ template void sparse_basic(const SparseMatrixType& re VERIFY_IS_APPROX(m1+=m2, refM1+=refM2); VERIFY_IS_APPROX(m1-=m2, refM1-=refM2); + VERIFY_IS_APPROX(m1.col(0).dot(refM2.row(0)), refM1.col(0).dot(refM2.row(0))); + refM4.setRandom(); // sparse cwise* dense VERIFY_IS_APPROX(m3.cwise()*refM4, refM3.cwise()*refM4); @@ -281,69 +283,6 @@ template void sparse_basic(const SparseMatrixType& re VERIFY_IS_APPROX(m2.transpose().eval(), refMat2.transpose().eval()); VERIFY_IS_APPROX(m2.transpose(), refMat2.transpose()); } - - // test matrix product - { - DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); - DenseMatrix refMat3 = DenseMatrix::Zero(rows, rows); - DenseMatrix refMat4 = DenseMatrix::Zero(rows, rows); - DenseMatrix dm4 = DenseMatrix::Zero(rows, rows); - SparseMatrixType m2(rows, rows); - SparseMatrixType m3(rows, rows); - SparseMatrixType m4(rows, rows); - initSparse(density, refMat2, m2); - initSparse(density, refMat3, m3); - initSparse(density, refMat4, m4); - VERIFY_IS_APPROX(m4=m2*m3, refMat4=refMat2*refMat3); - VERIFY_IS_APPROX(m4=m2.transpose()*m3, refMat4=refMat2.transpose()*refMat3); - VERIFY_IS_APPROX(m4=m2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose()); - VERIFY_IS_APPROX(m4=m2*m3.transpose(), refMat4=refMat2*refMat3.transpose()); - - // sparse * dense - VERIFY_IS_APPROX(dm4=m2*refMat3, refMat4=refMat2*refMat3); - VERIFY_IS_APPROX(dm4=m2*refMat3.transpose(), refMat4=refMat2*refMat3.transpose()); - VERIFY_IS_APPROX(dm4=m2.transpose()*refMat3, refMat4=refMat2.transpose()*refMat3); - VERIFY_IS_APPROX(dm4=m2.transpose()*refMat3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose()); - - // dense * sparse - VERIFY_IS_APPROX(dm4=refMat2*m3, refMat4=refMat2*refMat3); - VERIFY_IS_APPROX(dm4=refMat2*m3.transpose(), refMat4=refMat2*refMat3.transpose()); - VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3, refMat4=refMat2.transpose()*refMat3); - VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose()); - } - - // test self adjoint products - { - DenseMatrix b = DenseMatrix::Random(rows, rows); - DenseMatrix x = DenseMatrix::Random(rows, rows); - DenseMatrix refX = DenseMatrix::Random(rows, rows); - DenseMatrix refUp = DenseMatrix::Zero(rows, rows); - DenseMatrix refLo = DenseMatrix::Zero(rows, rows); - DenseMatrix refS = DenseMatrix::Zero(rows, rows); - SparseMatrixType mUp(rows, rows); - SparseMatrixType mLo(rows, rows); - SparseMatrixType mS(rows, rows); - do { - initSparse(density, refUp, mUp, ForceRealDiag|/*ForceNonZeroDiag|*/MakeUpperTriangular); - } while (refUp.isZero()); - refLo = refUp.transpose().conjugate(); - mLo = mUp.transpose().conjugate(); - refS = refUp + refLo; - refS.diagonal() *= 0.5; - mS = mUp + mLo; - for (int k=0; k()*b, refX=refS*b); - VERIFY_IS_APPROX(x=mLo.template marked()*b, refX=refS*b); - VERIFY_IS_APPROX(x=mS.template marked()*b, refX=refS*b); - } // test prune { diff --git a/test/sparse_product.cpp b/test/sparse_product.cpp new file mode 100644 index 000000000..95b3546c1 --- /dev/null +++ b/test/sparse_product.cpp @@ -0,0 +1,132 @@ +// 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 Daniel Gomez Ferro +// +// 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 . + +#include "sparse.h" + +template void sparse_product(const SparseMatrixType& ref) +{ + const int rows = ref.rows(); + const int cols = ref.cols(); + typedef typename SparseMatrixType::Scalar Scalar; + enum { Flags = SparseMatrixType::Flags }; + + double density = std::max(8./(rows*cols), 0.01); + typedef Matrix DenseMatrix; + typedef Matrix DenseVector; + Scalar eps = 1e-6; + + // test matrix-matrix product + /* + { + DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); + DenseMatrix refMat3 = DenseMatrix::Zero(rows, rows); + DenseMatrix refMat4 = DenseMatrix::Zero(rows, rows); + DenseMatrix dm4 = DenseMatrix::Zero(rows, rows); + SparseMatrixType m2(rows, rows); + SparseMatrixType m3(rows, rows); + SparseMatrixType m4(rows, rows); + initSparse(density, refMat2, m2); + initSparse(density, refMat3, m3); + initSparse(density, refMat4, m4); + VERIFY_IS_APPROX(m4=m2*m3, refMat4=refMat2*refMat3); + VERIFY_IS_APPROX(m4=m2.transpose()*m3, refMat4=refMat2.transpose()*refMat3); + VERIFY_IS_APPROX(m4=m2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose()); + VERIFY_IS_APPROX(m4=m2*m3.transpose(), refMat4=refMat2*refMat3.transpose()); + + // sparse * dense + VERIFY_IS_APPROX(dm4=m2*refMat3, refMat4=refMat2*refMat3); + VERIFY_IS_APPROX(dm4=m2*refMat3.transpose(), refMat4=refMat2*refMat3.transpose()); + VERIFY_IS_APPROX(dm4=m2.transpose()*refMat3, refMat4=refMat2.transpose()*refMat3); + VERIFY_IS_APPROX(dm4=m2.transpose()*refMat3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose()); + + // dense * sparse + VERIFY_IS_APPROX(dm4=refMat2*m3, refMat4=refMat2*refMat3); + VERIFY_IS_APPROX(dm4=refMat2*m3.transpose(), refMat4=refMat2*refMat3.transpose()); + VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3, refMat4=refMat2.transpose()*refMat3); + VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose()); + } + */ + + // test matrix - diagonal product + { + DenseMatrix refM2 = DenseMatrix::Zero(rows, rows); + DenseMatrix refM3 = DenseMatrix::Zero(rows, rows); + DiagonalMatrix d1(DenseVector::Random(rows)); + SparseMatrixType m2(rows, rows); + SparseMatrixType m3(rows, rows); + initSparse(density, refM2, m2); + initSparse(density, refM3, m3); +// std::cerr << "foo\n" << (m2*d1).toDense() << "\n\n" << refM2*d1 << "\n\n"; + VERIFY_IS_APPROX(m3=m2*d1, refM3=refM2*d1); + VERIFY_IS_APPROX(m3=m2.transpose()*d1, refM3=refM2.transpose()*d1); + VERIFY_IS_APPROX(m3=d1*m2, refM3=d1*refM2); +// std::cerr << "foo\n" << (d1*m2.transpose()).toDense() << "\n\n" << d1 * refM2.transpose() << "\n\n"; + VERIFY_IS_APPROX(m3=d1*m2.transpose(), refM3=d1 * refM2.transpose()); + } + + // test self adjoint products +// { +// DenseMatrix b = DenseMatrix::Random(rows, rows); +// DenseMatrix x = DenseMatrix::Random(rows, rows); +// DenseMatrix refX = DenseMatrix::Random(rows, rows); +// DenseMatrix refUp = DenseMatrix::Zero(rows, rows); +// DenseMatrix refLo = DenseMatrix::Zero(rows, rows); +// DenseMatrix refS = DenseMatrix::Zero(rows, rows); +// SparseMatrixType mUp(rows, rows); +// SparseMatrixType mLo(rows, rows); +// SparseMatrixType mS(rows, rows); +// do { +// initSparse(density, refUp, mUp, ForceRealDiag|/*ForceNonZeroDiag|*/MakeUpperTriangular); +// } while (refUp.isZero()); +// refLo = refUp.transpose().conjugate(); +// mLo = mUp.transpose().conjugate(); +// refS = refUp + refLo; +// refS.diagonal() *= 0.5; +// mS = mUp + mLo; +// for (int k=0; k()*b, refX=refS*b); +// VERIFY_IS_APPROX(x=mLo.template marked()*b, refX=refS*b); +// VERIFY_IS_APPROX(x=mS.template marked()*b, refX=refS*b); +// } + +} + +void test_sparse_product() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST( sparse_product(SparseMatrix(8, 8)) ); + CALL_SUBTEST( sparse_product(SparseMatrix >(16, 16)) ); + CALL_SUBTEST( sparse_product(SparseMatrix(33, 33)) ); + + CALL_SUBTEST( sparse_product(DynamicSparseMatrix(8, 8)) ); + } +} diff --git a/test/sparse_vector.cpp b/test/sparse_vector.cpp index 8207e522a..934719f2c 100644 --- a/test/sparse_vector.cpp +++ b/test/sparse_vector.cpp @@ -84,6 +84,7 @@ template void sparse_vector(int rows, int cols) VERIFY_IS_APPROX(v1-=v2, refV1-=refV2); VERIFY_IS_APPROX(v1.dot(v2), refV1.dot(refV2)); + VERIFY_IS_APPROX(v1.dot(refV2), refV1.dot(refV2)); }