| // This file is part of Eigen, a lightweight C++ template library |
| // for linear algebra. |
| // |
| // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> |
| // |
| // This Source Code Form is subject to the terms of the Mozilla |
| // Public License v. 2.0. If a copy of the MPL was not distributed |
| // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
| |
| #ifndef EIGEN_RANDOMSETTER_H |
| #define EIGEN_RANDOMSETTER_H |
| |
| #if defined(EIGEN_GOOGLEHASH_SUPPORT) |
| // Ensure the ::google namespace exists, required for checking existence of |
| // ::google::dense_hash_map and ::google::sparse_hash_map. |
| namespace google {} |
| #endif |
| |
| #include "./InternalHeaderCheck.h" |
| |
| namespace Eigen { |
| |
| /** Represents a std::map |
| * |
| * \see RandomSetter |
| */ |
| template<typename Scalar> struct StdMapTraits |
| { |
| typedef int KeyType; |
| typedef std::map<KeyType,Scalar> Type; |
| enum { |
| IsSorted = 1 |
| }; |
| |
| static void setInvalidKey(Type&, const KeyType&) {} |
| }; |
| |
| |
| /** Represents a std::unordered_map |
| * \see RandomSetter |
| */ |
| template<typename Scalar> struct StdUnorderedMapTraits |
| { |
| typedef int KeyType; |
| typedef std::unordered_map<KeyType,Scalar> Type; |
| enum { |
| IsSorted = 0 |
| }; |
| |
| static void setInvalidKey(Type&, const KeyType&) {} |
| }; |
| |
| #if defined(EIGEN_GOOGLEHASH_SUPPORT) |
| |
| namespace google { |
| |
| // Namespace work-around, since sometimes dense_hash_map and sparse_hash_map |
| // are in the global namespace, and other times they are under ::google. |
| using namespace ::google; |
| |
| template<typename KeyType, typename Scalar> |
| struct DenseHashMap { |
| typedef dense_hash_map<KeyType, Scalar> type; |
| }; |
| |
| template<typename KeyType, typename Scalar> |
| struct SparseHashMap { |
| typedef sparse_hash_map<KeyType, Scalar> type; |
| }; |
| |
| } // namespace google |
| |
| /** Represents a google::dense_hash_map |
| * |
| * \see RandomSetter |
| */ |
| template<typename Scalar> struct GoogleDenseHashMapTraits |
| { |
| typedef int KeyType; |
| typedef typename google::DenseHashMap<KeyType,Scalar>::type Type; |
| enum { |
| IsSorted = 0 |
| }; |
| |
| static void setInvalidKey(Type& map, const KeyType& k) |
| { map.set_empty_key(k); } |
| }; |
| |
| /** Represents a google::sparse_hash_map |
| * |
| * \see RandomSetter |
| */ |
| template<typename Scalar> struct GoogleSparseHashMapTraits |
| { |
| typedef int KeyType; |
| typedef typename google::SparseHashMap<KeyType,Scalar>::type Type; |
| enum { |
| IsSorted = 0 |
| }; |
| |
| static void setInvalidKey(Type&, const KeyType&) {} |
| }; |
| #endif |
| |
| /** \class RandomSetter |
| * \ingroup SparseExtra_Module |
| * \brief The RandomSetter is a wrapper object allowing to set/update a sparse matrix with random access |
| * |
| * \tparam SparseMatrixType the type of the sparse matrix we are updating |
| * \tparam MapTraits a traits class representing the map implementation used for the temporary sparse storage. |
| * Its default value depends on the system. |
| * \tparam OuterPacketBits defines the number of rows (or columns) manage by a single map object |
| * as a power of two exponent. |
| * |
| * This class temporarily represents a sparse matrix object using a generic map implementation allowing for |
| * efficient random access. The conversion from the compressed representation to a hash_map object is performed |
| * in the RandomSetter constructor, while the sparse matrix is updated back at destruction time. This strategy |
| * suggest the use of nested blocks as in this example: |
| * |
| * \code |
| * SparseMatrix<double> m(rows,cols); |
| * { |
| * RandomSetter<SparseMatrix<double> > w(m); |
| * // don't use m but w instead with read/write random access to the coefficients: |
| * for(;;) |
| * w(rand(),rand()) = rand; |
| * } |
| * // when w is deleted, the data are copied back to m |
| * // and m is ready to use. |
| * \endcode |
| * |
| * Since hash_map objects are not fully sorted, representing a full matrix as a single hash_map would |
| * involve a big and costly sort to update the compressed matrix back. To overcome this issue, a RandomSetter |
| * use multiple hash_map, each representing 2^OuterPacketBits columns or rows according to the storage order. |
| * To reach optimal performance, this value should be adjusted according to the average number of nonzeros |
| * per rows/columns. |
| * |
| * The possible values for the template parameter MapTraits are: |
| * - \b StdMapTraits: corresponds to std::map. (does not perform very well) |
| * - \b StdUnorderedMapTraits: corresponds to std::unordered_map |
| * - \b GoogleDenseHashMapTraits: corresponds to google::dense_hash_map (best efficiency, reasonable memory consumption) |
| * - \b GoogleSparseHashMapTraits: corresponds to google::sparse_hash_map (best memory consumption, relatively good performance) |
| * |
| * The default map implementation depends on the availability, and the preferred order is: |
| * GoogleSparseHashMapTraits, StdUnorderedMapTraits, and finally StdMapTraits. |
| * |
| * For performance and memory consumption reasons it is highly recommended to use one of |
| * Google's hash_map implementations. To enable the support for them, you must define |
| * EIGEN_GOOGLEHASH_SUPPORT. This will include both <google/dense_hash_map> and |
| * <google/sparse_hash_map> for you. |
| * |
| * \see https://github.com/sparsehash/sparsehash |
| */ |
| template<typename SparseMatrixType, |
| template <typename T> class MapTraits = |
| #if defined(EIGEN_GOOGLEHASH_SUPPORT) |
| GoogleDenseHashMapTraits |
| #else |
| StdUnorderedMapTraits |
| #endif |
| ,int OuterPacketBits = 6> |
| class RandomSetter |
| { |
| typedef typename SparseMatrixType::Scalar Scalar; |
| typedef typename SparseMatrixType::StorageIndex StorageIndex; |
| |
| struct ScalarWrapper |
| { |
| ScalarWrapper() : value(0) {} |
| Scalar value; |
| }; |
| typedef typename MapTraits<ScalarWrapper>::KeyType KeyType; |
| typedef typename MapTraits<ScalarWrapper>::Type HashMapType; |
| static const int OuterPacketMask = (1 << OuterPacketBits) - 1; |
| enum { |
| SwapStorage = 1 - MapTraits<ScalarWrapper>::IsSorted, |
| TargetRowMajor = (SparseMatrixType::Flags & RowMajorBit) ? 1 : 0, |
| SetterRowMajor = SwapStorage ? 1-TargetRowMajor : TargetRowMajor |
| }; |
| |
| public: |
| |
| /** Constructs a random setter object from the sparse matrix \a target |
| * |
| * Note that the initial value of \a target are imported. If you want to re-set |
| * a sparse matrix from scratch, then you must set it to zero first using the |
| * setZero() function. |
| */ |
| inline RandomSetter(SparseMatrixType& target) |
| : mp_target(&target) |
| { |
| const Index outerSize = SwapStorage ? target.innerSize() : target.outerSize(); |
| const Index innerSize = SwapStorage ? target.outerSize() : target.innerSize(); |
| m_outerPackets = outerSize >> OuterPacketBits; |
| if (outerSize&OuterPacketMask) |
| m_outerPackets += 1; |
| m_hashmaps = new HashMapType[m_outerPackets]; |
| // compute number of bits needed to store inner indices |
| Index aux = innerSize - 1; |
| m_keyBitsOffset = 0; |
| while (aux) |
| { |
| ++m_keyBitsOffset; |
| aux = aux >> 1; |
| } |
| KeyType ik = (1<<(OuterPacketBits+m_keyBitsOffset)); |
| for (Index k=0; k<m_outerPackets; ++k) |
| MapTraits<ScalarWrapper>::setInvalidKey(m_hashmaps[k],ik); |
| |
| // insert current coeffs |
| for (Index j=0; j<mp_target->outerSize(); ++j) |
| for (typename SparseMatrixType::InnerIterator it(*mp_target,j); it; ++it) |
| (*this)(TargetRowMajor?j:it.index(), TargetRowMajor?it.index():j) = it.value(); |
| } |
| |
| /** Destructor updating back the sparse matrix target */ |
| ~RandomSetter() |
| { |
| KeyType keyBitsMask = (1<<m_keyBitsOffset)-1; |
| if (!SwapStorage) // also means the map is sorted |
| { |
| mp_target->setZero(); |
| mp_target->makeCompressed(); |
| mp_target->reserve(nonZeros()); |
| Index prevOuter = -1; |
| for (Index k=0; k<m_outerPackets; ++k) |
| { |
| const Index outerOffset = (1<<OuterPacketBits) * k; |
| typename HashMapType::iterator end = m_hashmaps[k].end(); |
| for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it!=end; ++it) |
| { |
| const Index outer = (it->first >> m_keyBitsOffset) + outerOffset; |
| const Index inner = it->first & keyBitsMask; |
| if (prevOuter!=outer) |
| { |
| for (Index j=prevOuter+1;j<=outer;++j) |
| mp_target->startVec(j); |
| prevOuter = outer; |
| } |
| mp_target->insertBackByOuterInner(outer, inner) = it->second.value; |
| } |
| } |
| mp_target->finalize(); |
| } |
| else |
| { |
| VectorXi positions(mp_target->outerSize()); |
| positions.setZero(); |
| // pass 1 |
| for (Index k=0; k<m_outerPackets; ++k) |
| { |
| typename HashMapType::iterator end = m_hashmaps[k].end(); |
| for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it!=end; ++it) |
| { |
| const Index outer = it->first & keyBitsMask; |
| ++positions[outer]; |
| } |
| } |
| // prefix sum |
| StorageIndex count = 0; |
| for (Index j=0; j<mp_target->outerSize(); ++j) |
| { |
| StorageIndex tmp = positions[j]; |
| mp_target->outerIndexPtr()[j] = count; |
| positions[j] = count; |
| count += tmp; |
| } |
| mp_target->makeCompressed(); |
| mp_target->outerIndexPtr()[mp_target->outerSize()] = count; |
| mp_target->resizeNonZeros(count); |
| // pass 2 |
| for (Index k=0; k<m_outerPackets; ++k) |
| { |
| const Index outerOffset = (1<<OuterPacketBits) * k; |
| typename HashMapType::iterator end = m_hashmaps[k].end(); |
| for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it!=end; ++it) |
| { |
| const Index inner = (it->first >> m_keyBitsOffset) + outerOffset; |
| const Index outer = it->first & keyBitsMask; |
| // sorted insertion |
| // Note that we have to deal with at most 2^OuterPacketBits unsorted coefficients, |
| // moreover those 2^OuterPacketBits coeffs are likely to be sparse, an so only a |
| // small fraction of them have to be sorted, whence the following simple procedure: |
| Index posStart = mp_target->outerIndexPtr()[outer]; |
| Index i = (positions[outer]++) - 1; |
| while ( (i >= posStart) && (mp_target->innerIndexPtr()[i] > inner) ) |
| { |
| mp_target->valuePtr()[i+1] = mp_target->valuePtr()[i]; |
| mp_target->innerIndexPtr()[i+1] = mp_target->innerIndexPtr()[i]; |
| --i; |
| } |
| mp_target->innerIndexPtr()[i+1] = internal::convert_index<StorageIndex>(inner); |
| mp_target->valuePtr()[i+1] = it->second.value; |
| } |
| } |
| } |
| delete[] m_hashmaps; |
| } |
| |
| /** \returns a reference to the coefficient at given coordinates \a row, \a col */ |
| Scalar& operator() (Index row, Index col) |
| { |
| const Index outer = SetterRowMajor ? row : col; |
| const Index inner = SetterRowMajor ? col : row; |
| const Index outerMajor = outer >> OuterPacketBits; // index of the packet/map |
| const Index outerMinor = outer & OuterPacketMask; // index of the inner vector in the packet |
| const KeyType key = internal::convert_index<KeyType>((outerMinor<<m_keyBitsOffset) | inner); |
| return m_hashmaps[outerMajor][key].value; |
| } |
| |
| /** \returns the number of non zero coefficients |
| * |
| * \note According to the underlying map/hash_map implementation, |
| * this function might be quite expensive. |
| */ |
| Index nonZeros() const |
| { |
| Index nz = 0; |
| for (Index k=0; k<m_outerPackets; ++k) |
| nz += static_cast<Index>(m_hashmaps[k].size()); |
| return nz; |
| } |
| |
| |
| protected: |
| |
| HashMapType* m_hashmaps; |
| SparseMatrixType* mp_target; |
| Index m_outerPackets; |
| unsigned char m_keyBitsOffset; |
| }; |
| |
| } // end namespace Eigen |
| |
| #endif // EIGEN_RANDOMSETTER_H |
