From e5d301dc961ddfaba6e38c497904b2aee378a7cc Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 26 Jun 2008 23:22:26 +0000 Subject: [PATCH] various work on the Sparse module: * added some glue to Eigen/Core (SparseBit, ei_eval, Matrix) * add two new sparse matrix types: HashMatrix: based on std::map (for random writes) LinkedVectorMatrix: array of linked vectors (for outer coherent writes, e.g. to transpose a matrix) * add a SparseSetter class to easily set/update any kind of matrices, e.g.: { SparseSetter wrapper(mymatrix); for (...) wrapper->coeffRef(rand(),rand()) = rand(); } * automatic shallow copy for RValue * and a lot of mess ! plus: * remove the remaining ArrayBit related stuff * don't use alloca in product for very large memory allocation --- Eigen/Sparse | 8 + Eigen/src/Core/CacheFriendlyProduct.h | 10 +- Eigen/src/Core/Matrix.h | 3 +- Eigen/src/Core/Product.h | 1 - Eigen/src/Core/Transpose.h | 2 + Eigen/src/Core/util/Constants.h | 19 +- Eigen/src/Core/util/Meta.h | 4 +- Eigen/src/Sparse/CoreIterators.h | 41 ++-- Eigen/src/Sparse/HashMatrix.h | 165 +++++++++++++++ Eigen/src/Sparse/LinkedVectorMatrix.h | 288 ++++++++++++++++++++++++++ Eigen/src/Sparse/SparseArray.h | 32 ++- Eigen/src/Sparse/SparseMatrix.h | 198 ++++++------------ Eigen/src/Sparse/SparseMatrixBase.h | 163 +++++++++++++++ Eigen/src/Sparse/SparseProduct.h | 169 +++++++++++++++ Eigen/src/Sparse/SparseSetter.h | 102 +++++++++ Eigen/src/Sparse/SparseUtil.h | 61 ++++++ bench/sparse_01.cpp | 129 +++++++++++- 17 files changed, 1226 insertions(+), 169 deletions(-) create mode 100644 Eigen/src/Sparse/HashMatrix.h create mode 100644 Eigen/src/Sparse/LinkedVectorMatrix.h create mode 100644 Eigen/src/Sparse/SparseMatrixBase.h create mode 100644 Eigen/src/Sparse/SparseProduct.h create mode 100644 Eigen/src/Sparse/SparseSetter.h create mode 100644 Eigen/src/Sparse/SparseUtil.h diff --git a/Eigen/Sparse b/Eigen/Sparse index 6a5eaf6a0..962660cbb 100644 --- a/Eigen/Sparse +++ b/Eigen/Sparse @@ -4,12 +4,20 @@ #include "Core" #include #include +#include +#include namespace Eigen { +#include "src/Sparse/SparseUtil.h" +#include "src/Sparse/SparseMatrixBase.h" +#include "src/Sparse/SparseProduct.h" #include "src/Sparse/SparseArray.h" #include "src/Sparse/SparseMatrix.h" +#include "src/Sparse/HashMatrix.h" +#include "src/Sparse/LinkedVectorMatrix.h" #include "src/Sparse/CoreIterators.h" +#include "src/Sparse/SparseSetter.h" } // namespace Eigen diff --git a/Eigen/src/Core/CacheFriendlyProduct.h b/Eigen/src/Core/CacheFriendlyProduct.h index 4a0e4e24a..0da84eeac 100644 --- a/Eigen/src/Core/CacheFriendlyProduct.h +++ b/Eigen/src/Core/CacheFriendlyProduct.h @@ -86,7 +86,12 @@ static void ei_cache_friendly_product( const int l2BlockRows = MaxL2BlockSize > rows ? rows : MaxL2BlockSize; const int l2BlockCols = MaxL2BlockSize > cols ? cols : MaxL2BlockSize; const int l2BlockSize = MaxL2BlockSize > size ? size : MaxL2BlockSize; - Scalar* __restrict__ block = (Scalar*)alloca(sizeof(Scalar)*l2BlockRows*size); + Scalar* __restrict__ block = 0; + const int allocBlockSize = sizeof(Scalar)*l2BlockRows*size; + if (allocBlockSize>16000000) + block = (Scalar*)malloc(allocBlockSize); + else + block = (Scalar*)alloca(allocBlockSize); Scalar* __restrict__ rhsCopy = (Scalar*)alloca(sizeof(Scalar)*l2BlockSize); // loops on each L2 cache friendly blocks of the result @@ -347,6 +352,9 @@ static void ei_cache_friendly_product( } } } + + if (allocBlockSize>16000000) + free(block); } #endif // EIGEN_CACHE_FRIENDLY_PRODUCT_H diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h index 40947e997..f00a27b33 100644 --- a/Eigen/src/Core/Matrix.h +++ b/Eigen/src/Core/Matrix.h @@ -92,7 +92,8 @@ struct ei_traits > _Rows, _Cols, _MaxRows, _MaxCols, _Flags >::ret, - CoeffReadCost = NumTraits::ReadCost + CoeffReadCost = NumTraits::ReadCost, + SupportedAccessPatterns = RandomAccessPattern }; }; diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index a0b29449c..e3b4b18e2 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -258,7 +258,6 @@ template inline const typename ProductReturnType::Type MatrixBase::operator*(const MatrixBase &other) const { - assert( (Derived::Flags&ArrayBit) == (OtherDerived::Flags&ArrayBit) ); return typename ProductReturnType::Type(derived(), other.derived()); } diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h index f19d4affa..c364e5f86 100644 --- a/Eigen/src/Core/Transpose.h +++ b/Eigen/src/Core/Transpose.h @@ -63,6 +63,8 @@ template class Transpose EIGEN_GENERIC_PUBLIC_INTERFACE(Transpose) + class InnerIterator; + inline Transpose(const MatrixType& matrix) : m_matrix(matrix) {} EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Transpose) diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 716d86243..eafcbca7f 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -138,10 +138,8 @@ const unsigned int LowerTriangularBit = 0x400; /** \ingroup flags * - * means the object is just an array of scalars, and operations on it are regarded as operations - * on every of these scalars taken separately. - */ -const unsigned int ArrayBit = 0x800; + * means the expression includes sparse matrices and the sparse path has to be taken. */ +const unsigned int SparseBit = 0x800; /** \ingroup flags * @@ -155,7 +153,7 @@ const unsigned int HereditaryBits = RowMajorBit | EvalBeforeNestingBit | EvalBeforeAssigningBit | LargeBit - | ArrayBit; + | SparseBit; // Possible values for the Mode parameter of part() and of extract() const unsigned int Upper = UpperTriangularBit; @@ -173,7 +171,7 @@ enum { Aligned=0, UnAligned=1 }; enum { ConditionalJumpCost = 5 }; enum CornerType { TopLeft, TopRight, BottomLeft, BottomRight }; enum DirectionType { Vertical, Horizontal }; -enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, DiagonalProduct }; +enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, DiagonalProduct, SparseProduct }; enum { InnerVectorization, @@ -188,5 +186,14 @@ enum { NoUnrolling }; +enum { + Dense = 0, + Sparse = SparseBit +}; + +const int FullyCoherentAccessPattern = 0x1; +const int InnerCoherentAccessPattern = 0x2 | FullyCoherentAccessPattern; +const int OuterCoherentAccessPattern = 0x4 | InnerCoherentAccessPattern; +const int RandomAccessPattern = 0x8 | OuterCoherentAccessPattern; #endif // EIGEN_CONSTANTS_H diff --git a/Eigen/src/Core/util/Meta.h b/Eigen/src/Core/util/Meta.h index 5d809f619..93d50441e 100644 --- a/Eigen/src/Core/util/Meta.h +++ b/Eigen/src/Core/util/Meta.h @@ -175,7 +175,9 @@ template struct ei_size_at_compile_time enum { ret = (_Rows==Dynamic || _Cols==Dynamic) ? Dynamic : _Rows * _Cols }; }; -template class ei_eval +template::Flags&SparseBit> class ei_eval; + +template class ei_eval { typedef typename ei_traits::Scalar _Scalar; enum {_Rows = ei_traits::RowsAtCompileTime, diff --git a/Eigen/src/Sparse/CoreIterators.h b/Eigen/src/Sparse/CoreIterators.h index 00f3ae781..f44db0c1f 100644 --- a/Eigen/src/Sparse/CoreIterators.h +++ b/Eigen/src/Sparse/CoreIterators.h @@ -25,30 +25,47 @@ #ifndef EIGEN_COREITERATORS_H #define EIGEN_COREITERATORS_H +/* This file contains the respective InnerIterator definition of the expressions defined in Eigen/Core + */ + template class MatrixBase::InnerIterator { typedef typename Derived::Scalar Scalar; public: - InnerIterator(const Derived& mat, int col) - : m_matrix(mat), m_row(0), m_col(col), m_end(mat.rows()) + InnerIterator(const Derived& mat, int outer) + : m_matrix(mat), m_inner(0), m_outer(outer), m_end(mat.rows()) {} - Scalar value() { return m_matrix.coeff(m_row, m_col); } + Scalar value() + { + return (Derived::Flags&RowMajorBit) ? m_matrix.coeff(m_outer, m_inner) + : m_matrix.coeff(m_inner, m_outer); + } - InnerIterator& operator++() { m_row++; return *this; } + InnerIterator& operator++() { m_inner++; return *this; } - int index() const { return m_row; } + int index() const { return m_inner; } - operator bool() const { return m_row < m_end && m_row>=0; } + operator bool() const { return m_inner < m_end && m_inner>=0; } protected: const Derived& m_matrix; - int m_row; - const int m_col; + int m_inner; + const int m_outer; const int m_end; }; +template +class Transpose::InnerIterator : public MatrixType::InnerIterator +{ + public: + + InnerIterator(const Transpose& trans, int outer) + : MatrixType::InnerIterator(trans.m_matrix, outer) + {} +}; + template class CwiseUnaryOp::InnerIterator { @@ -57,8 +74,8 @@ class CwiseUnaryOp::InnerIterator typedef typename _MatrixTypeNested::InnerIterator MatrixTypeIterator; public: - InnerIterator(const CwiseUnaryOp& unaryOp, int col) - : m_iter(unaryOp.m_matrix,col), m_functor(unaryOp.m_functor), m_id(-1) + InnerIterator(const CwiseUnaryOp& unaryOp, int outer) + : m_iter(unaryOp.m_matrix,outer), m_functor(unaryOp.m_functor), m_id(-1) { this->operator++(); } @@ -101,8 +118,8 @@ class CwiseBinaryOp::InnerIterator typedef typename _RhsNested::InnerIterator RhsIterator; public: - InnerIterator(const CwiseBinaryOp& binOp, int col) - : m_lhsIter(binOp.m_lhs,col), m_rhsIter(binOp.m_rhs,col), m_functor(binOp.m_functor), m_id(-1) + InnerIterator(const CwiseBinaryOp& binOp, int outer) + : m_lhsIter(binOp.m_lhs,outer), m_rhsIter(binOp.m_rhs,outer), m_functor(binOp.m_functor), m_id(-1) { this->operator++(); } diff --git a/Eigen/src/Sparse/HashMatrix.h b/Eigen/src/Sparse/HashMatrix.h new file mode 100644 index 000000000..753c7c5e9 --- /dev/null +++ b/Eigen/src/Sparse/HashMatrix.h @@ -0,0 +1,165 @@ +// 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_HASHMATRIX_H +#define EIGEN_HASHMATRIX_H + +template +struct ei_traits > +{ + typedef _Scalar Scalar; + enum { + RowsAtCompileTime = Dynamic, + ColsAtCompileTime = Dynamic, + MaxRowsAtCompileTime = Dynamic, + MaxColsAtCompileTime = Dynamic, + Flags = _Flags, + CoeffReadCost = NumTraits::ReadCost, + SupportedAccessPatterns = RandomAccessPattern + }; +}; + +// TODO reimplement this class using custom linked lists +template +class HashMatrix : public SparseMatrixBase > +{ + public: + EIGEN_GENERIC_PUBLIC_INTERFACE(HashMatrix) + class InnerIterator; + protected: + + typedef typename std::map::iterator MapIterator; + typedef typename std::map::const_iterator ConstMapIterator; + + public: + inline int rows() const { return m_innerSize; } + inline int cols() const { return m_data.size(); } + + inline const Scalar& coeff(int row, int col) const + { + const MapIterator it = m_data[col].find(row); + if (it!=m_data[col].end()) + return Scalar(0); + return it->second; + } + + inline Scalar& coeffRef(int row, int col) + { + return m_data[col][row]; + } + + public: + + inline void startFill(int reserveSize = 1000) {} + + inline Scalar& fill(int row, int col) { return coeffRef(row, col); } + + inline void endFill() {} + + ~HashMatrix() + {} + + inline void shallowCopy(const HashMatrix& other) + { + EIGEN_DBG_SPARSE(std::cout << "HashMatrix:: shallowCopy\n"); + // FIXME implement a true shallow copy !! + resize(other.rows(), other.cols()); + for (int j=0; jouterSize(); ++j) + m_data[j] = other.m_data[j]; + } + + void resize(int _rows, int _cols) + { + if (cols() != _cols) + { + m_data.resize(_cols); + } + m_innerSize = _rows; + } + + inline HashMatrix(int rows, int cols) + : m_innerSize(0) + { + resize(rows, cols); + } + + template + inline HashMatrix(const MatrixBase& other) + : m_innerSize(0) + { + *this = other.derived(); + } + + inline HashMatrix& operator=(const HashMatrix& other) + { + if (other.isRValue()) + { + shallowCopy(other); + } + else + { + resize(other.rows(), other.cols()); + for (int col=0; col + inline HashMatrix& operator=(const MatrixBase& other) + { + return SparseMatrixBase::operator=(other); + } + + protected: + + std::vector > m_data; + int m_innerSize; + +}; + +template +class HashMatrix::InnerIterator +{ + public: + + InnerIterator(const HashMatrix& mat, int col) + : m_matrix(mat), m_it(mat.m_data[col].begin()), m_end(mat.m_data[col].end()) + {} + + InnerIterator& operator++() { m_it++; return *this; } + + Scalar value() { return m_it->second; } + + int index() const { return m_it->first; } + + operator bool() const { return m_it!=m_end; } + + protected: + const HashMatrix& m_matrix; + typename HashMatrix::ConstMapIterator m_it; + typename HashMatrix::ConstMapIterator m_end; +}; + +#endif // EIGEN_HASHMATRIX_H diff --git a/Eigen/src/Sparse/LinkedVectorMatrix.h b/Eigen/src/Sparse/LinkedVectorMatrix.h new file mode 100644 index 000000000..2e2618e26 --- /dev/null +++ b/Eigen/src/Sparse/LinkedVectorMatrix.h @@ -0,0 +1,288 @@ +// 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_LINKEDVECTORMATRIX_H +#define EIGEN_LINKEDVECTORMATRIX_H + +template +struct ei_traits > +{ + typedef _Scalar Scalar; + enum { + RowsAtCompileTime = Dynamic, + ColsAtCompileTime = Dynamic, + MaxRowsAtCompileTime = Dynamic, + MaxColsAtCompileTime = Dynamic, + Flags = _Flags, + CoeffReadCost = NumTraits::ReadCost, + SupportedAccessPatterns = InnerCoherentAccessPattern + }; +}; + +template +struct LinkedVector +{ + LinkedVector() : size(0), next(0), prev(0) {} + Element data[BlockSize]; + LinkedVector* next; + LinkedVector* prev; + int size; + bool isFull() const { return size==BlockSize; } +}; + +template +class LinkedVectorMatrix : public SparseMatrixBase > +{ + public: + EIGEN_GENERIC_PUBLIC_INTERFACE(LinkedVectorMatrix) + class InnerIterator; + protected: + + enum { + RowMajor = Flags&RowMajorBit ? 1 : 0 + }; + + struct ValueIndex + { + ValueIndex() : value(0), index(0) {} + ValueIndex(Scalar v, int i) : value(v), index(i) {} + Scalar value; + int index; + }; + typedef LinkedVector LinkedVectorBlock; + + inline int find(LinkedVectorBlock** _el, int id) + { + LinkedVectorBlock* el = *_el; + while (el && el->data[el->size-1].indexnext; + *_el = el; + if (el) + { + // binary search + int maxI = el->size-1; + int minI = 0; + int i = el->size/2; + const ValueIndex* data = el->data; + while (data[i].index!=id) + { + if (data[i].index=maxI) + return -1; + } + if (data[i].index==id) + return i; + } + return -1; + } + + public: + inline int rows() const { return RowMajor ? m_data.size() : m_innerSize; } + inline int cols() const { return RowMajor ? m_innerSize : m_data.size(); } + + inline const Scalar& coeff(int row, int col) const + { + const int outer = RowMajor ? row : col; + const int inner = RowMajor ? col : row; + + LinkedVectorBlock* el = m_data[outer]; + int id = find(&el, inner); + if (id<0) + return Scalar(0); + return el->data[id].value; + } + + inline Scalar& coeffRef(int row, int col) + { + const int outer = RowMajor ? row : col; + const int inner = RowMajor ? col : row; + + LinkedVectorBlock* el = m_data[outer]; + int id = find(&el, inner); + ei_assert(id>=0); +// if (id<0) +// return Scalar(0); + return el->data[id].value; + } + + public: + + inline void startFill(int reserveSize = 1000) + { + clear(); + for (int i=0; idata[m_ends[outer]->size-1].index < inner); + if (m_ends[outer]->isFull()) + { + + LinkedVectorBlock* el = new LinkedVectorBlock(); + m_ends[outer]->next = el; + el->prev = m_ends[outer]; + m_ends[outer] = el; + } + } + m_ends[outer]->data[m_ends[outer]->size].index = inner; + return m_ends[outer]->data[m_ends[outer]->size++].value; + } + + inline void endFill() + { + } + + ~LinkedVectorMatrix() + { + if (this->isNotShared()) + clear(); + } + + void clear() + { + for (int i=0; inext; + delete tmp; + } + } + } + + void resize(int rows, int cols) + { + const int outers = RowMajor ? rows : cols; + const int inners = RowMajor ? cols : rows; + + if (this->outerSize() != outers) + { + clear(); + m_data.resize(outers); + m_ends.resize(outers); + for (int i=0; i + inline LinkedVectorMatrix(const MatrixBase& other) + : m_innerSize(0) + { + *this = other.derived(); + } + + inline void shallowCopy(const LinkedVectorMatrix& other) + { + EIGEN_DBG_SPARSE(std::cout << "LinkedVectorMatrix:: shallowCopy\n"); + resize(other.rows(), other.cols()); + for (int j=0; jouterSize(); ++j) + { + m_data[j] = other.m_data[j]; + m_ends[j] = other.m_ends[j]; + } + other.markAsCopied(); + } + + inline LinkedVectorMatrix& operator=(const LinkedVectorMatrix& other) + { + if (other.isRValue()) + { + shallowCopy(other); + return *this; + } + else + { + // TODO implement a specialized deep copy here + return operator=(other); + } + } + + template + inline LinkedVectorMatrix& operator=(const MatrixBase& other) + { + return SparseMatrixBase::operator=(other.derived()); + } + + protected: + + std::vector m_data; + std::vector m_ends; + int m_innerSize; + +}; + + +template +class LinkedVectorMatrix::InnerIterator +{ + public: + + InnerIterator(const LinkedVectorMatrix& mat, int col) + : m_matrix(mat), m_el(mat.m_data[col]), m_it(0) + {} + + InnerIterator& operator++() { if (m_itsize) m_it++; else {m_el = m_el->next; m_it=0;}; return *this; } + + Scalar value() { return m_el->data[m_it].value; } + + int index() const { return m_el->data[m_it].index; } + + operator bool() const { return m_el && (m_el->next || m_itsize); } + + protected: + const LinkedVectorMatrix& m_matrix; + LinkedVectorBlock* m_el; + int m_it; +}; + +#endif // EIGEN_LINKEDVECTORMATRIX_H diff --git a/Eigen/src/Sparse/SparseArray.h b/Eigen/src/Sparse/SparseArray.h index 2a87801a8..3ac181e04 100644 --- a/Eigen/src/Sparse/SparseArray.h +++ b/Eigen/src/Sparse/SparseArray.h @@ -32,11 +32,11 @@ template class SparseArray { public: SparseArray() - : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0) + : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0), m_isShared(0) {} SparseArray(int size) - : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0) + : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0), m_isShared(0) { resize(size); } @@ -51,6 +51,28 @@ template class SparseArray resize(other.size()); memcpy(m_values, other.m_values, m_size * sizeof(Scalar)); memcpy(m_indices, other.m_indices, m_size * sizeof(int)); + m_isShared = 0; + } + + void shallowCopy(const SparseArray& other) + { + delete[] m_values; + delete[] m_indices; + m_values = other.m_values; + m_indices = other.m_indices; + m_size = other.m_size; + m_allocatedSize = other.m_allocatedSize; + m_isShared = false; + other.m_isShared = true; + } + + ~SparseArray() + { + if (!m_isShared) + { + delete[] m_values; + delete[] m_indices; + } } void reserve(int size) @@ -117,7 +139,11 @@ template class SparseArray Scalar* m_values; int* m_indices; int m_size; - int m_allocatedSize; + struct { + int m_allocatedSize:31; + mutable int m_isShared:1; + }; + }; #endif // EIGEN_SPARSE_ARRAY_H diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h index d3ee6e3b1..d9d404b8a 100644 --- a/Eigen/src/Sparse/SparseMatrix.h +++ b/Eigen/src/Sparse/SparseMatrix.h @@ -25,8 +25,6 @@ #ifndef EIGEN_SPARSEMATRIX_H #define EIGEN_SPARSEMATRIX_H -template class SparseMatrix; - /** \class SparseMatrix * * \brief Sparse matrix @@ -36,8 +34,8 @@ template class SparseMatrix; * See http://www.netlib.org/linalg/html_templates/node91.html for details on the storage scheme. * */ -template -struct ei_traits > +template +struct ei_traits > { typedef _Scalar Scalar; enum { @@ -45,13 +43,16 @@ struct ei_traits > ColsAtCompileTime = Dynamic, MaxRowsAtCompileTime = Dynamic, MaxColsAtCompileTime = Dynamic, - Flags = 0, - CoeffReadCost = NumTraits::ReadCost + Flags = _Flags, + CoeffReadCost = NumTraits::ReadCost, + SupportedAccessPatterns = FullyCoherentAccessPattern }; }; -template -class SparseMatrix : public MatrixBase > + + +template +class SparseMatrix : public SparseMatrixBase > { public: EIGEN_GENERIC_PUBLIC_INTERFACE(SparseMatrix) @@ -63,6 +64,8 @@ class SparseMatrix : public MatrixBase > int m_rows; int m_cols; + public: + inline int rows() const { return m_rows; } inline int cols() const { return m_cols; } @@ -81,8 +84,8 @@ class SparseMatrix : public MatrixBase > inline Scalar& coeffRef(int row, int col) { - int id = m_colPtrs[cols]; - int end = m_colPtrs[cols+1]; + int id = m_colPtrs[col]; + int end = m_colPtrs[col+1]; while (id > class InnerIterator; - inline int rows() const { return rows(); } - inline int cols() const { return cols(); } /** \returns the number of non zero coefficients */ inline int nonZeros() const { return m_data.size(); } - inline const Scalar& operator() (int row, int col) const - { - return coeff(row, col); - } - - inline Scalar& operator() (int row, int col) - { - return coeffRef(row, col); - } - inline void startFill(int reserveSize = 1000) { m_data.clear(); @@ -170,143 +161,80 @@ class SparseMatrix : public MatrixBase > resize(rows, cols); } + inline void shallowCopy(const SparseMatrix& other) + { + EIGEN_DBG_SPARSE(std::cout << "SparseMatrix:: shallowCopy\n"); + delete[] m_colPtrs; + m_colPtrs = 0; + m_rows = other.rows(); + m_cols = other.cols(); + m_colPtrs = other.m_colPtrs; + m_data.shallowCopy(other.m_data); + other.markAsCopied(); + } + inline SparseMatrix& operator=(const SparseMatrix& other) { - resize(other.rows(), other.cols()); - m_colPtrs = other.m_colPtrs; - for (int col=0; col<=cols(); ++col) - m_colPtrs[col] = other.m_colPtrs[col]; - m_data = other.m_data; - return *this; + if (other.isRValue()) + { + shallowCopy(other); + } + else + { + resize(other.rows(), other.cols()); + for (int col=0; col<=cols(); ++col) + m_colPtrs[col] = other.m_colPtrs[col]; + m_data = other.m_data; + return *this; + } } template inline SparseMatrix& operator=(const MatrixBase& other) { - resize(other.rows(), other.cols()); - startFill(std::max(m_rows,m_cols)*2); - for (int col=0; col::operator=(other); } + template + SparseMatrix operator*(const MatrixBase& other) + { + SparseMatrix res(rows(), other.cols()); + ei_sparse_product(*this,other.derived(),res); + return res; + } - // old explicit operator+ -// template -// SparseMatrix operator+(const Other& other) -// { -// SparseMatrix res(rows(), cols()); -// res.startFill(nonZeros()*3); -// for (int col=0; col()); -// } - - - // WARNING for efficiency reason it currently outputs the transposed matrix friend std::ostream & operator << (std::ostream & s, const SparseMatrix& m) { - s << "Nonzero entries:\n"; - for (uint i=0; i&>(m); return s; } /** Destructor */ inline ~SparseMatrix() { - delete[] m_colPtrs; + if (this->isNotShared()) + delete[] m_colPtrs; } }; -template -class SparseMatrix::InnerIterator +template +class SparseMatrix::InnerIterator { public: InnerIterator(const SparseMatrix& mat, int col) diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h new file mode 100644 index 000000000..6357cbeff --- /dev/null +++ b/Eigen/src/Sparse/SparseMatrixBase.h @@ -0,0 +1,163 @@ +// 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_SPARSEMATRIXBASE_H +#define EIGEN_SPARSEMATRIXBASE_H + +template +class SparseMatrixBase : public MatrixBase +{ + public: + + typedef MatrixBase Base; + typedef typename Base::Scalar Scalar; + enum { + Flags = Base::Flags, + RowMajor = ei_traits::Flags&RowMajorBit ? 1 : 0 + }; + + 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)); } + + SparseMatrixBase() + : m_isRValue(false), m_hasBeenCopied(false) + {} + + bool isRValue() const { return m_isRValue; } + Derived& temporary() { m_isRValue = true; return derived(); } + + int outerSize() const { return RowMajor ? this->rows() : this->cols(); } + int innerSize() const { return RowMajor ? this->cols() : this->rows(); } + + inline Derived& operator=(const Derived& other) + { + if (other.isRValue()) + { + m_hasBeenCopied = true; + derived().shallowCopy(other); + } + else + this->operator=(other); + return derived(); + } + + + + template + inline Derived& operator=(const MatrixBase& other) + { + const bool transpose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit); + const int outerSize = (int(OtherDerived::Flags) & RowMajorBit) ? other.rows() : other.cols(); + // eval to a temporary and then do a shallow copy + typename ei_meta_if, Derived>::ret temp(other.rows(), other.cols()); + + temp.startFill(std::max(this->rows(),this->cols())*2); +// std::cout << other.rows() << " xm " << other.cols() << "\n"; + for (int j=0; j + inline Derived& operator=(const SparseMatrixBase& other) + { + const bool transpose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit); + const int outerSize = (int(OtherDerived::Flags) & RowMajorBit) ? other.rows() : other.cols(); + if ((!transpose) && other.isRValue()) + { + // eval without temporary + derived().resize(other.rows(), other.cols()); + derived().startFill(std::max(this->rows(),this->cols())*2); + for (int j=0; joperator=(static_cast&>(other)); + return derived(); + } + + friend std::ostream & operator << (std::ostream & s, const SparseMatrixBase& m) + { + + if (Flags&RowMajorBit) + { + for (int row=0; row trans = m.derived(); + s << trans; + } + return s; + } + + protected: + + bool isNotShared() { return !m_hasBeenCopied; } + void markAsCopied() const { m_hasBeenCopied = true; } + + protected: + + bool m_isRValue; + mutable bool m_hasBeenCopied; +}; + + + +#endif // EIGEN_SPARSEMATRIXBASE_H diff --git a/Eigen/src/Sparse/SparseProduct.h b/Eigen/src/Sparse/SparseProduct.h new file mode 100644 index 000000000..9abddbd27 --- /dev/null +++ b/Eigen/src/Sparse/SparseProduct.h @@ -0,0 +1,169 @@ +// 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_SPARSEPRODUCT_H +#define EIGEN_SPARSEPRODUCT_H + +#define DENSE_TMP 1 +#define MAP_TMP 2 +#define LIST_TMP 3 + +#define TMP_TMP 3 + +template +struct ListEl +{ + int next; + int index; + Scalar value; +}; + +template +static void ei_sparse_product(const Lhs& lhs, const Rhs& rhs, SparseMatrix::Scalar>& res) +{ + int rows = lhs.rows(); + int cols = rhs.rows(); + int size = lhs.cols(); + + float ratio = std::max(float(lhs.nonZeros())/float(lhs.rows()*lhs.cols()), float(rhs.nonZeros())/float(rhs.rows()*rhs.cols())); + std::cout << ratio << "\n"; + + ei_assert(size == rhs.rows()); + typedef typename ei_traits::Scalar Scalar; + #if (TMP_TMP == MAP_TMP) + std::map tmp; + #elif (TMP_TMP == LIST_TMP) + std::vector > tmp(2*rows); + #else + std::vector tmp(rows); + #endif + res.resize(rows, cols); + res.startFill(2*std::max(rows, cols)); + for (int j=0; j::iterator hint = tmp.begin(); + typename std::map::iterator r; + #elif (TMP_TMP == LIST_TMP) + int tmp_el = tmp_start; + #endif + for (typename Lhs::InnerIterator lhsIt(lhs, rhsIt.index()); lhsIt; ++lhsIt) + { + #if (TMP_TMP == MAP_TMP) + r = hint; + Scalar v = lhsIt.value() * rhsIt.value(); + int id = lhsIt.index(); + while (r!=tmp.end() && r->first < id) + ++r; + if (r!=tmp.end() && r->first==id) + { + r->second += v; + hint = r; + } + else + hint = tmp.insert(r, std::pair(id, v)); + ++hint; + #elif (TMP_TMP == LIST_TMP) + Scalar v = lhsIt.value() * rhsIt.value(); + int id = lhsIt.index(); + if (tmp_size==0) + { + tmp_start = 0; + tmp_el = 0; + tmp_size++; + tmp[0].value = v; + tmp[0].index = id; + tmp[0].next = -1; + } + else if (id= 0 && tmp[nextel].index<=id) + { + tmp_el = nextel; + nextel = tmp[nextel].next; + } + + if (tmp[tmp_el].index==id) + { + tmp[tmp_el].value += v; + } + else + { + tmp[tmp_size].value = v; + tmp[tmp_size].index = id; + tmp[tmp_size].next = tmp[tmp_el].next; + tmp[tmp_el].next = tmp_size; + tmp_size++; + } + } + #else + tmp[lhsIt.index()] += lhsIt.value() * rhsIt.value(); + #endif + //res.coeffRef(lhsIt.index(), j) += lhsIt.value() * rhsIt.value(); + } + } + #if (TMP_TMP == MAP_TMP) + for (typename std::map::const_iterator k=tmp.begin(); k!=tmp.end(); ++k) + if (k->second!=0) + res.fill(k->first, j) = k->second; + #elif (TMP_TMP == LIST_TMP) + int k = tmp_start; + while (k>=0) + { + if (tmp[k].value!=0) + res.fill(tmp[k].index, j) = tmp[k].value; + k = tmp[k].next; + } + #else + for (int k=0; k " << float(res.nonZeros())/float(res.rows()*res.cols()) << "\n"; +} + +#endif // EIGEN_SPARSEPRODUCT_H diff --git a/Eigen/src/Sparse/SparseSetter.h b/Eigen/src/Sparse/SparseSetter.h new file mode 100644 index 000000000..93b634322 --- /dev/null +++ b/Eigen/src/Sparse/SparseSetter.h @@ -0,0 +1,102 @@ +// 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_SPARSESETTER_H +#define EIGEN_SPARSESETTER_H + +template::ret> +struct ei_sparse_setter_selector; + +template::type> +class SparseSetter +{ + typedef typename ei_unref::type _WrapperType; + public: + + inline SparseSetter(MatrixType& matrix) : m_wrapper(matrix), mp_matrix(&matrix) {} + + ~SparseSetter() + { *mp_matrix = m_wrapper; } + + inline _WrapperType* operator->() { return &m_wrapper; } + + inline _WrapperType& operator*() { return m_wrapper; } + + protected: + + WrapperType m_wrapper; + MatrixType* mp_matrix; +}; + +template +struct ei_sparse_setter_selector +{ + typedef MatrixType& type; +}; + +// forward each derived of SparseMatrixBase to the generic SparseMatrixBase specializations +template +struct ei_sparse_setter_selector, AccessPattern, AccessPatternNotSupported> +: public ei_sparse_setter_selector >,AccessPattern, AccessPatternNotSupported> +{}; + +template +struct ei_sparse_setter_selector, AccessPattern, AccessPatternNotSupported> +: public ei_sparse_setter_selector >,AccessPattern, AccessPatternNotSupported> +{}; + +template +struct ei_sparse_setter_selector, AccessPattern, AccessPatternNotSupported> +: public ei_sparse_setter_selector >,AccessPattern, AccessPatternNotSupported> +{}; + +// generic SparseMatrixBase specializations +template +struct ei_sparse_setter_selector, RandomAccessPattern, AccessPatternNotSupported> +{ + typedef HashMatrix type; +}; + +template +struct ei_sparse_setter_selector, OuterCoherentAccessPattern, AccessPatternNotSupported> +{ + typedef HashMatrix type; +}; + +template +struct ei_sparse_setter_selector, InnerCoherentAccessPattern, AccessPatternNotSupported> +{ + typedef LinkedVectorMatrix type; +}; + +template +struct ei_sparse_setter_selector, FullyCoherentAccessPattern, AccessPatternNotSupported> +{ + typedef SparseMatrix type; +}; + +#endif // EIGEN_SPARSESETTER_H diff --git a/Eigen/src/Sparse/SparseUtil.h b/Eigen/src/Sparse/SparseUtil.h new file mode 100644 index 000000000..e89a7c322 --- /dev/null +++ b/Eigen/src/Sparse/SparseUtil.h @@ -0,0 +1,61 @@ +// 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_SPARSEUTIL_H +#define EIGEN_SPARSEUTIL_H + +#ifdef NDEBUG +#define EIGEN_DBG_SPARSE(X) +#else +#define EIGEN_DBG_SPARSE(X) X +#endif + +template class SparseMatrix; +template class HashMatrix; +template class LinkedVectorMatrix; + +const int AccessPatternNotSupported = 0x0; +const int AccessPatternSupported = 0x1; + + +template struct ei_support_access_pattern +{ + enum { ret = (int(ei_traits::SupportedAccessPatterns) & AccessPattern) == AccessPattern + ? AccessPatternSupported + : AccessPatternNotSupported + }; +}; + +template class ei_eval +{ + typedef typename ei_traits::Scalar _Scalar; + enum { + _Flags = ei_traits::Flags + }; + + public: + typedef SparseMatrix<_Scalar, _Flags> type; +}; + +#endif // EIGEN_SPARSEUTIL_H diff --git a/bench/sparse_01.cpp b/bench/sparse_01.cpp index 9521a17cf..69ab4dda3 100644 --- a/bench/sparse_01.cpp +++ b/bench/sparse_01.cpp @@ -12,15 +12,23 @@ using namespace Eigen; USING_PART_OF_NAMESPACE_EIGEN #ifndef REPEAT -#define REPEAT 40000000 +#define REPEAT 10 #endif +#define REPEATPRODUCT 1 + +#define SIZE 10 +#define DENSITY 0.2 + +// #define NODENSEMATRIX + typedef MatrixXf DenseMatrix; +// typedef Matrix DenseMatrix; typedef SparseMatrix EigenSparseMatrix; typedef gmm::csc_matrix GmmSparse; typedef gmm::col_matrix< gmm::wsvector > GmmDynSparse; -void fillMatrix(float density, int rows, int cols, MatrixXf* pDenseMatrix, EigenSparseMatrix* pSparseMatrix, GmmSparse* pGmmMatrix=0) +void fillMatrix(float density, int rows, int cols, DenseMatrix* pDenseMatrix, EigenSparseMatrix* pSparseMatrix, GmmSparse* pGmmMatrix=0) { GmmDynSparse gmmT(rows, cols); if (pSparseMatrix) @@ -49,52 +57,155 @@ void fillMatrix(float density, int rows, int cols, MatrixXf* pDenseMatrix, Eigen int main(int argc, char *argv[]) { - int rows = 4000; - int cols = 4000; - float density = 0.1; + int rows = SIZE; + int cols = SIZE; + float density = DENSITY; // dense matrices + #ifndef NODENSEMATRIX DenseMatrix m1(rows,cols), m2(rows,cols), m3(rows,cols), m4(rows,cols); + #endif // sparse matrices EigenSparseMatrix sm1(rows,cols), sm2(rows,cols), sm3(rows,cols), sm4(rows,cols); + HashMatrix hm4(rows,cols); + // GMM++ matrices GmmDynSparse gmmT4(rows,cols); GmmSparse gmmM1(rows,cols), gmmM2(rows,cols), gmmM3(rows,cols), gmmM4(rows,cols); + #ifndef NODENSEMATRIX fillMatrix(density, rows, cols, &m1, &sm1, &gmmM1); fillMatrix(density, rows, cols, &m2, &sm2, &gmmM2); fillMatrix(density, rows, cols, &m3, &sm3, &gmmM3); + #else + fillMatrix(density, rows, cols, 0, &sm1, &gmmM1); + fillMatrix(density, rows, cols, 0, &sm2, &gmmM2); + fillMatrix(density, rows, cols, 0, &sm3, &gmmM3); + #endif BenchTimer timer; + //-------------------------------------------------------------------------------- + // COEFF WISE OPERATORS + //-------------------------------------------------------------------------------- +#if 1 + std::cout << "\n\n\"m4 = m1 + m2 + 2 * m3\":\n\n"; + + timer.reset(); timer.start(); - for (int k=0; k<10; ++k) + asm("#begin"); + for (int k=0; k lm5(rows, cols); + lm5 = lm4; + lm5 = sm4; + cout << endl << lm5 << endl; + + sm3 = sm4.transpose(); + cout << endl << lm5 << endl; + + cout << endl << "SM1 before random editing: " << endl << sm1 << endl; + { + SparseSetter w1(sm1); + w1->coeffRef(4,2) = ei_random(); + w1->coeffRef(2,6) = ei_random(); + w1->coeffRef(0,4) = ei_random(); + w1->coeffRef(9,3) = ei_random(); + } + cout << endl << "SM1 after random editing: " << endl << sm1 << endl; +#endif + return 0; }