From ea1990ef3d95a2e042b0ececdc4f21c0f5473cc2 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 23 Jun 2008 13:25:22 +0000 Subject: [PATCH] add experimental code for sparse matrix: - uses the common "Compressed Column Storage" scheme - supports every unary and binary operators with xpr template assuming binaryOp(0,0) == 0 and unaryOp(0) = 0 (otherwise a sparse matrix doesnot make sense) - this is the first commit, so of course, there are still several shorcommings ! --- Eigen/Sparse | 16 ++ Eigen/src/Core/CwiseBinaryOp.h | 2 + Eigen/src/Core/CwiseUnaryOp.h | 2 + Eigen/src/Core/MatrixBase.h | 2 + Eigen/src/Sparse/CoreIterators.h | 152 ++++++++++++++ Eigen/src/Sparse/SparseArray.h | 123 ++++++++++++ Eigen/src/Sparse/SparseMatrix.h | 331 +++++++++++++++++++++++++++++++ bench/sparse_01.cpp | 100 ++++++++++ 8 files changed, 728 insertions(+) create mode 100644 Eigen/Sparse create mode 100644 Eigen/src/Sparse/CoreIterators.h create mode 100644 Eigen/src/Sparse/SparseArray.h create mode 100644 Eigen/src/Sparse/SparseMatrix.h create mode 100644 bench/sparse_01.cpp diff --git a/Eigen/Sparse b/Eigen/Sparse new file mode 100644 index 000000000..6a5eaf6a0 --- /dev/null +++ b/Eigen/Sparse @@ -0,0 +1,16 @@ +#ifndef EIGEN_SPARSE_MODULE_H +#define EIGEN_SPARSE_MODULE_H + +#include "Core" +#include +#include + +namespace Eigen { + +#include "src/Sparse/SparseArray.h" +#include "src/Sparse/SparseMatrix.h" +#include "src/Sparse/CoreIterators.h" + +} // namespace Eigen + +#endif // EIGEN_SPARSE_MODULE_H diff --git a/Eigen/src/Core/CwiseBinaryOp.h b/Eigen/src/Core/CwiseBinaryOp.h index b12aa65bd..ec4619781 100644 --- a/Eigen/src/Core/CwiseBinaryOp.h +++ b/Eigen/src/Core/CwiseBinaryOp.h @@ -84,6 +84,8 @@ class CwiseBinaryOp : ei_no_assignment_operator, typedef typename ei_traits::LhsNested LhsNested; typedef typename ei_traits::RhsNested RhsNested; + class InnerIterator; + inline CwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp()) : m_lhs(lhs), m_rhs(rhs), m_functor(func) { diff --git a/Eigen/src/Core/CwiseUnaryOp.h b/Eigen/src/Core/CwiseUnaryOp.h index 10228e3ec..7c466d8c2 100644 --- a/Eigen/src/Core/CwiseUnaryOp.h +++ b/Eigen/src/Core/CwiseUnaryOp.h @@ -69,6 +69,8 @@ class CwiseUnaryOp : ei_no_assignment_operator, EIGEN_GENERIC_PUBLIC_INTERFACE(CwiseUnaryOp) + class InnerIterator; + inline CwiseUnaryOp(const MatrixType& mat, const UnaryOp& func = UnaryOp()) : m_matrix(mat), m_functor(func) {} diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index ad0f0b1e8..a31f2028f 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -55,6 +55,8 @@ template class MatrixBase public: + class InnerIterator; + typedef typename ei_traits::Scalar Scalar; typedef typename ei_packet_traits::type PacketScalar; diff --git a/Eigen/src/Sparse/CoreIterators.h b/Eigen/src/Sparse/CoreIterators.h new file mode 100644 index 000000000..00f3ae781 --- /dev/null +++ b/Eigen/src/Sparse/CoreIterators.h @@ -0,0 +1,152 @@ +// 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_COREITERATORS_H +#define EIGEN_COREITERATORS_H + +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()) + {} + + Scalar value() { return m_matrix.coeff(m_row, m_col); } + + InnerIterator& operator++() { m_row++; return *this; } + + int index() const { return m_row; } + + operator bool() const { return m_row < m_end && m_row>=0; } + + protected: + const Derived& m_matrix; + int m_row; + const int m_col; + const int m_end; +}; + +template +class CwiseUnaryOp::InnerIterator +{ + typedef typename CwiseUnaryOp::Scalar Scalar; + typedef typename ei_traits::_MatrixTypeNested _MatrixTypeNested; + 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) + { + this->operator++(); + } + + InnerIterator& operator++() + { + if (m_iter) + { + m_id = m_iter.index(); + m_value = m_functor(m_iter.value()); + ++m_iter; + } + else + { + m_id = -1; + } + return *this; + } + + Scalar value() const { return m_value; } + + int index() const { return m_id; } + + operator bool() const { return m_id>=0; } + + protected: + MatrixTypeIterator m_iter; + const UnaryOp& m_functor; + Scalar m_value; + int m_id; +}; + +template +class CwiseBinaryOp::InnerIterator +{ + typedef typename CwiseBinaryOp::Scalar Scalar; + typedef typename ei_traits::_LhsNested _LhsNested; + typedef typename _LhsNested::InnerIterator LhsIterator; + typedef typename ei_traits::_RhsNested _RhsNested; + typedef typename _RhsNested::InnerIterator RhsIterator; + public: + + 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) + { + this->operator++(); + } + + InnerIterator& operator++() + { + if (m_lhsIter && m_rhsIter && (m_lhsIter.index() == m_rhsIter.index())) + { + m_id = m_lhsIter.index(); + m_value = m_functor(m_lhsIter.value(), m_rhsIter.value()); + ++m_lhsIter; + ++m_rhsIter; + } + else if (m_lhsIter && ((!m_rhsIter) || m_lhsIter.index() < m_rhsIter.index())) + { + m_id = m_lhsIter.index(); + m_value = m_functor(m_lhsIter.value(), Scalar(0)); + ++m_lhsIter; + } + else if (m_rhsIter && ((!m_lhsIter) || m_lhsIter.index() > m_rhsIter.index())) + { + m_id = m_rhsIter.index(); + m_value = m_functor(Scalar(0), m_rhsIter.value()); + ++m_rhsIter; + } + else + { + m_id = -1; + } + return *this; + } + + Scalar value() const { return m_value; } + + int index() const { return m_id; } + + operator bool() const { return m_id>=0; } + + protected: + LhsIterator m_lhsIter; + RhsIterator m_rhsIter; + const BinaryOp& m_functor; + Scalar m_value; + int m_id; +}; + +#endif // EIGEN_COREITERATORS_H diff --git a/Eigen/src/Sparse/SparseArray.h b/Eigen/src/Sparse/SparseArray.h new file mode 100644 index 000000000..2a87801a8 --- /dev/null +++ b/Eigen/src/Sparse/SparseArray.h @@ -0,0 +1,123 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2008 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_SPARSE_ARRAY_H +#define EIGEN_SPARSE_ARRAY_H + +/** Stores a sparse set of values as a list of values and a list of indices. + * + */ +template class SparseArray +{ + public: + SparseArray() + : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0) + {} + + SparseArray(int size) + : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0) + { + resize(size); + } + + SparseArray(const SparseArray& other) + { + *this = other; + } + + SparseArray& operator=(const SparseArray& other) + { + resize(other.size()); + memcpy(m_values, other.m_values, m_size * sizeof(Scalar)); + memcpy(m_indices, other.m_indices, m_size * sizeof(int)); + } + + void reserve(int size) + { + int newAllocatedSize = m_size + size; + if (newAllocatedSize > m_allocatedSize) + { + Scalar* newValues = new Scalar[newAllocatedSize]; + int* newIndices = new int[newAllocatedSize]; + // copy + memcpy(newValues, m_values, m_size * sizeof(Scalar)); + memcpy(newIndices, m_indices, m_size * sizeof(int)); + // delete old stuff + delete[] m_values; + delete[] m_indices; + m_values = newValues; + m_indices = newIndices; + m_allocatedSize = newAllocatedSize; + } + } + + void resize(int size, int reserveSizeFactor = 0) + { + if (m_allocatedSize +// +// 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_SPARSEMATRIX_H +#define EIGEN_SPARSEMATRIX_H + +template class SparseMatrix; + +/** \class SparseMatrix + * + * \brief Sparse matrix + * + * \param _Scalar the scalar type, i.e. the type of the coefficients + * + * See http://www.netlib.org/linalg/html_templates/node91.html for details on the storage scheme. + * + */ +template +struct ei_traits > +{ + typedef _Scalar Scalar; + enum { + RowsAtCompileTime = Dynamic, + ColsAtCompileTime = Dynamic, + MaxRowsAtCompileTime = Dynamic, + MaxColsAtCompileTime = Dynamic, + Flags = 0, + CoeffReadCost = NumTraits::ReadCost + }; +}; + +template +class SparseMatrix : public MatrixBase > +{ + public: + EIGEN_GENERIC_PUBLIC_INTERFACE(SparseMatrix) + + protected: + + int* m_colPtrs; + SparseArray m_data; + int m_rows; + int m_cols; + + inline int _rows() const { return m_rows; } + inline int _cols() const { return m_cols; } + + inline const Scalar& _coeff(int row, int col) const + { + int id = m_colPtrs[col]; + int end = m_colPtrs[col+1]; + while (id=0 && m_colPtrs[i]==0) + { + m_colPtrs[i] = m_data.size(); + --i; + } + m_colPtrs[col+1] = m_colPtrs[col]; + } + assert(m_colPtrs[col+1] == m_data.size()); + int id = m_colPtrs[col+1]; + m_colPtrs[col+1]++; + + m_data.append(0, row); + return m_data.value(id); + } + + inline void endFill() + { + int size = m_data.size(); + int i = m_cols; + // find the last filled column + while (i>=0 && m_colPtrs[i]==0) + --i; + i++; + while (i<=m_cols) + { + m_colPtrs[i] = size; + ++i; + } + } + + void resize(int rows, int cols) + { + if (m_cols != cols) + { + delete[] m_colPtrs; + m_colPtrs = new int [cols+1]; + m_rows = rows; + m_cols = cols; + } + } + + inline SparseMatrix(int rows, int cols) + : m_rows(0), m_cols(0), m_colPtrs(0) + { + resize(rows, cols); + } + + 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; + } + + 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 +// 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 +class SparseMatrix::InnerIterator +{ + public: + InnerIterator(const SparseMatrix& mat, int col) + : m_matrix(mat), m_id(mat.m_colPtrs[col]), m_start(m_id), m_end(mat.m_colPtrs[col+1]) + {} + + InnerIterator& operator++() { m_id++; return *this; } + + Scalar value() { return m_matrix.m_data.value(m_id); } + + int index() const { return m_matrix.m_data.index(m_id); } + + operator bool() const { return (m_id < m_end) && (m_id>=m_start); } + + protected: + const SparseMatrix& m_matrix; + int m_id; + const int m_start; + const int m_end; +}; + +#endif // EIGEN_SPARSEMATRIX_H diff --git a/bench/sparse_01.cpp b/bench/sparse_01.cpp new file mode 100644 index 000000000..9521a17cf --- /dev/null +++ b/bench/sparse_01.cpp @@ -0,0 +1,100 @@ + +// g++ -O3 -DNDEBUG sparse_01.cpp -I .. -o sparse_01 && ./sparse_01 + +#include +#include +#include + +#include "gmm/gmm.h" + +using namespace std; +using namespace Eigen; +USING_PART_OF_NAMESPACE_EIGEN + +#ifndef REPEAT +#define REPEAT 40000000 +#endif + +typedef MatrixXf 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) +{ + GmmDynSparse gmmT(rows, cols); + if (pSparseMatrix) + pSparseMatrix->startFill(rows*cols*density); + for(int j = 0; j < cols; j++) + { + for(int i = 0; i < rows; i++) + { + float v = (ei_random(0,1) < density) ? ei_random() : 0; + if (pDenseMatrix) + (*pDenseMatrix)(i,j) = v; + if (v!=0) + { + if (pSparseMatrix) + pSparseMatrix->fill(i,j) = v; + if (pGmmMatrix) + gmmT(i,j) = v; + } + } + } + if (pSparseMatrix) + pSparseMatrix->endFill(); + if (pGmmMatrix) + gmm::copy(gmmT, *pGmmMatrix); +} + +int main(int argc, char *argv[]) +{ + int rows = 4000; + int cols = 4000; + float density = 0.1; + + // dense matrices + DenseMatrix m1(rows,cols), m2(rows,cols), m3(rows,cols), m4(rows,cols); + + // sparse matrices + EigenSparseMatrix sm1(rows,cols), sm2(rows,cols), sm3(rows,cols), sm4(rows,cols); + + // GMM++ matrices + + GmmDynSparse gmmT4(rows,cols); + GmmSparse gmmM1(rows,cols), gmmM2(rows,cols), gmmM3(rows,cols), gmmM4(rows,cols); + + fillMatrix(density, rows, cols, &m1, &sm1, &gmmM1); + fillMatrix(density, rows, cols, &m2, &sm2, &gmmM2); + fillMatrix(density, rows, cols, &m3, &sm3, &gmmM3); + + BenchTimer timer; + + timer.start(); + for (int k=0; k<10; ++k) + m4 = m1 + m2 + 2 * m3; + timer.stop(); + std::cout << "Eigen dense = " << timer.value() << endl; + + timer.reset(); + timer.start(); + for (int k=0; k<10; ++k) + sm4 = sm1 + sm2 + 2 * sm3; + timer.stop(); + std::cout << "Eigen sparse = " << timer.value() << endl; + + timer.reset(); + timer.start(); + for (int k=0; k<10; ++k) + { + gmm::add(gmmM1, gmmM2, gmmT4); + gmm::add(gmm::scaled(gmmM3,2), gmmT4); + } + timer.stop(); + std::cout << "GMM++ sparse = " << timer.value() << endl; + +// sm3 = sm1 + m2; +// cout << m4.transpose() << "\n\n" << sm4 << endl; + return 0; +} +