#ifndef _utl_SparseMatrix_h_ #define _utl_SparseMatrix_h_ #include #include #include namespace utl { /** \class SparseMatrix SparseMatrix.h "utl/SparseMatrix.h" \brief Sparse container class for matrix data Importatant notes:
  • the operator[i][j] access is not creating entries. It should be used for reading of the whole matrix.
  • the operator(i, j) access creates the corresponding entry on the fly and initializes it to GetInitializer() before the reference to the element is returned
Please refer to the testSparseMatrix.cc and test SparseMatrixVectorOps.cc files for detailed usage cases. \author Darko Veberic \version $Id: SparseMatrix.h 23590 2013-05-06 17:59:58Z paul $ \date 5 Jun 2008 \ingroup math */ template class SparseMatrix { public: typedef int IndexType; /// row proxy class Row { public: Row(const SparseMatrix& m, const IndexType row) : fSparseMatrix(m), fRow(row) { } const T& operator[](const IndexType column) const { return fSparseMatrix(fRow, column); } private: const SparseMatrix& fSparseMatrix; const IndexType fRow; }; /// index class class Index { public: Index(const IndexType row, const IndexType column) : fRow(row), fColumn(column) { } IndexType GetRow() const { return fRow; } IndexType GetColumn() const { return fColumn; } bool operator<(const Index& i) const { return fRow < i.fRow || (fRow == i.fRow && fColumn < i.fColumn); } bool operator==(const Index& i) const { return fRow == i.fRow && fColumn == i.fColumn; } private: const IndexType fRow; const IndexType fColumn; }; private: typedef typename SparseMatrix::Index IndexKeyType; typedef std::map MapType; typedef typename MapType::value_type MapValueType; typedef typename MapType::iterator MapIterator; typedef typename MapType::const_iterator MapConstIterator; public: /// iterator class template class IteratorTransformer { public: IteratorTransformer(const IteratorType& it) : fIterator(it) { } IndexType GetRow() const { return fIterator->first.GetRow(); } IndexType GetColumn() const { return fIterator->first.GetColumn(); } ValueType& operator()() const { return fIterator->second; } ValueType& operator*() const { return (*this)(); } IteratorTransformer& operator++() { ++fIterator; return *this; } bool operator==(const IteratorTransformer& it) const { return fIterator == it.fIterator; } bool operator!=(const IteratorTransformer& it) const { return fIterator != it.fIterator; } private: IteratorType fIterator; }; SparseMatrix(const T& init = T()) : fInitializer(init) { } template explicit SparseMatrix(const SparseMatrix& m) : fInitializer(m.GetInitializer()) { for (typename SparseMatrix::ConstIterator it = m.SparseBegin(); it != m.SparseEnd(); ++it) (*this)(it.GetRow(), it.GetColumn()) = it(); } unsigned int GetNonEmptySize() const { return fMatrix.size(); } const T& GetInitializer() const { return fInitializer; } void Clear() { fMatrix.clear(); } /// noncreational [] access trough proxy typename SparseMatrix::Row operator[](const IndexType row) const { return typename SparseMatrix::Row(*this, row); } /// creational () access T& operator()(const IndexType row, const IndexType column) { return fMatrix.insert(MapValueType(IndexKeyType(row, column), fInitializer)).first->second; } /// noncreational () access for const const T& operator()(const IndexType row, const IndexType column) const { const MapConstIterator it = fMatrix.find(IndexKeyType(row, column)); return (it == fMatrix.end()) ? fInitializer : it->second; } bool IsEmpty() const { return fMatrix.empty(); } bool IsEmpty(const IndexType row, const IndexType column) const { const MapConstIterator it = fMatrix.find(IndexKeyType(row, column)); return it == fMatrix.end(); } unsigned int Reclaim() { unsigned int erased = 0; for (MapIterator it = fMatrix.begin(); it != fMatrix.end(); ) if (it->second != fInitializer) ++it; else { Erase(it++); ++erased; } return erased; } /// be sure to call Reclaim() first bool operator==(const SparseMatrix& m) const { return fMatrix == m.fMatrix && fInitializer == m.fInitializer; } /// be sure to call Reclaim() first bool operator!=(const SparseMatrix& m) const { return !operator==(m); } template bool operator==(const SparseMatrix& m) const { if (fInitializer != m.GetInitializer()) return false; ConstIterator it1 = SparseBegin(); typename SparseMatrix::ConstIterator it2 = m.SparseBegin(); while (it1 != SparseEnd() && it2 != m.SparseEnd()) { if (it1.GetRow() != it2.GetRow() || it1.GetColumn() != it2.GetColumn() || it1() != it2()) return false; ++it1; ++it2; } if (it1 != SparseEnd() || it2 != m.SparseEnd()) return false; return true; } template bool operator!=(const SparseMatrix& m) const { return !operator==(m); } typedef IteratorTransformer Iterator; typedef IteratorTransformer ConstIterator; Iterator SparseBegin() { return Iterator(fMatrix.begin()); } ConstIterator SparseBegin() const { return ConstIterator(fMatrix.begin()); } Iterator SparseEnd() { return Iterator(fMatrix.end()); } ConstIterator SparseEnd() const { return ConstIterator(fMatrix.end()); } Iterator Find(const IndexType row, const IndexType column) { return Iterator(fMatrix.find(Index(row, column))); } ConstIterator Find(const IndexType row, const IndexType column) const { return ConstIterator(fMatrix.find(Index(row, column))); } /// multiplication with scalar template SparseMatrix& operator*=(const U& value) { for (MapIterator it = fMatrix.begin(); it != fMatrix.end(); ++it) it->second *= value; return *this; } /// matrix addition template SparseMatrix& operator+=(const SparseMatrix& m) { for (typename SparseMatrix::ConstIterator it = m.SparseBegin(); it != m.SparseEnd(); ++it) (*this)(it.GetRow(), it.GetColumn()) += it(); return *this; } /// matrix subtraction template SparseMatrix& operator-=(const SparseMatrix& m) { for (typename SparseMatrix::ConstIterator it = m.SparseBegin(); it != m.SparseEnd(); ++it) (*this)(it.GetRow(), it.GetColumn()) -= it(); return *this; } /// in-situ transpose void Transpose() { MapType t; for (MapConstIterator it = fMatrix.begin(); it != fMatrix.end(); ++it) { const IndexKeyType& i = it->first; t.insert(t.begin(), std::make_pair(IndexKeyType(i.GetColumn(), i.GetRow()), it->second)); } fMatrix.swap(t); } private: void Erase(const MapIterator it) { fMatrix.erase(it); } MapType fMatrix; const T fInitializer; template friend class SparseMatrix; }; /** general matrix multiplication In case matrix elements are not simple scalar values, the corresponding operator* will be used. For example, imagine matrix elements are vectors and the element * element multiplication is equal to scalar product. The required dot product functor is then used. */ template typename boost::lambda::return_type_2< boost::lambda::arithmetic_action, SparseMatrix, SparseMatrix >::type operator*(const SparseMatrix& a, const SparseMatrix& b) { typedef typename boost::lambda::return_type_2< boost::lambda::arithmetic_action, SparseMatrix, SparseMatrix >::type ResultMatrix; ResultMatrix c; for (typename SparseMatrix::ConstIterator bIt = b.SparseBegin(); bIt != b.SparseEnd(); ++bIt) { const typename SparseMatrix::IndexType k = bIt.GetRow(); for (typename SparseMatrix::ConstIterator aIt = a.SparseBegin(); aIt != a.SparseEnd(); ++aIt) { const typename SparseMatrix::IndexType kk = aIt.GetColumn(); if (kk == k) { const typename SparseMatrix::IndexType i = aIt.GetRow(); const typename SparseMatrix::IndexType j = bIt.GetColumn(); c(i, j) += aIt() * bIt(); } } } return c; } // output stuff template Stream& operator<<(Stream& s, const SparseMatrix& m) { s << '('; typename SparseMatrix::ConstIterator it = m.SparseBegin(); if (it != m.SparseEnd()) { s << '{' << it.GetRow() << ',' << it.GetColumn() << ", " << it() << '}'; ++it; while (it != m.SparseEnd()) s << " {" << it.GetRow() << ',' << it.GetColumn() << ", " << it() << '}'; } return s << ')'; } // output functions assume the lowest index starts with 0, which is not // necessarily true template Stream& OutputSparse(Stream& s, const SparseMatrix& m, const typename SparseMatrix::IndexType rows, const typename SparseMatrix::IndexType columns, const char separator, const SparseType& sparse) { typedef typename SparseMatrix::IndexType Index; for (Index i = 0; i < rows; ++i) { if (m.IsEmpty(i, 0)) s << sparse; else s << m[i][0]; for (Index j = 1; j < columns; ++j) { s << separator; if (m.IsEmpty(i, j)) s << sparse; else s << m[i][j]; } s << '\n'; } return s; } template inline Stream& OutputSparse(Stream& s, const SparseMatrix& m, const typename SparseMatrix::IndexType rows, const typename SparseMatrix::IndexType columns) { return OutputSparse(s, m, rows, columns, ' ', '.'); } template inline Stream& Output(Stream& s, const SparseMatrix& m, const typename SparseMatrix::IndexType rows, const typename SparseMatrix::IndexType columns, const char separator = ' ') { return OutputSparse(s, m, rows, columns, separator, m.GetInitializer()); } } namespace boost { namespace lambda { /// specialization for custom class SMatrix<> template class plain_return_type_2< arithmetic_action, utl::SparseMatrix, utl::SparseMatrix > { private: // find type of T * U typedef typename return_type_2< arithmetic_action, T, U >::type res_type; public: typedef utl::SparseMatrix type; }; } } #include #endif