[136] | 1 | // This file is part of Eigen, a lightweight C++ template library
|
---|
| 2 | // for linear algebra.
|
---|
| 3 | //
|
---|
| 4 | // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
|
---|
| 5 | //
|
---|
| 6 | // This Source Code Form is subject to the terms of the Mozilla
|
---|
| 7 | // Public License v. 2.0. If a copy of the MPL was not distributed
|
---|
| 8 | // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
---|
| 9 |
|
---|
| 10 | #ifndef EIGEN_RANDOMSETTER_H
|
---|
| 11 | #define EIGEN_RANDOMSETTER_H
|
---|
| 12 |
|
---|
| 13 | namespace Eigen {
|
---|
| 14 |
|
---|
| 15 | /** Represents a std::map
|
---|
| 16 | *
|
---|
| 17 | * \see RandomSetter
|
---|
| 18 | */
|
---|
| 19 | template<typename Scalar> struct StdMapTraits
|
---|
| 20 | {
|
---|
| 21 | typedef int KeyType;
|
---|
| 22 | typedef std::map<KeyType,Scalar> Type;
|
---|
| 23 | enum {
|
---|
| 24 | IsSorted = 1
|
---|
| 25 | };
|
---|
| 26 |
|
---|
| 27 | static void setInvalidKey(Type&, const KeyType&) {}
|
---|
| 28 | };
|
---|
| 29 |
|
---|
| 30 | #ifdef EIGEN_UNORDERED_MAP_SUPPORT
|
---|
| 31 | /** Represents a std::unordered_map
|
---|
| 32 | *
|
---|
| 33 | * To use it you need to both define EIGEN_UNORDERED_MAP_SUPPORT and include the unordered_map header file
|
---|
| 34 | * yourself making sure that unordered_map is defined in the std namespace.
|
---|
| 35 | *
|
---|
| 36 | * For instance, with current version of gcc you can either enable C++0x standard (-std=c++0x) or do:
|
---|
| 37 | * \code
|
---|
| 38 | * #include <tr1/unordered_map>
|
---|
| 39 | * #define EIGEN_UNORDERED_MAP_SUPPORT
|
---|
| 40 | * namespace std {
|
---|
| 41 | * using std::tr1::unordered_map;
|
---|
| 42 | * }
|
---|
| 43 | * \endcode
|
---|
| 44 | *
|
---|
| 45 | * \see RandomSetter
|
---|
| 46 | */
|
---|
| 47 | template<typename Scalar> struct StdUnorderedMapTraits
|
---|
| 48 | {
|
---|
| 49 | typedef int KeyType;
|
---|
| 50 | typedef std::unordered_map<KeyType,Scalar> Type;
|
---|
| 51 | enum {
|
---|
| 52 | IsSorted = 0
|
---|
| 53 | };
|
---|
| 54 |
|
---|
| 55 | static void setInvalidKey(Type&, const KeyType&) {}
|
---|
| 56 | };
|
---|
| 57 | #endif // EIGEN_UNORDERED_MAP_SUPPORT
|
---|
| 58 |
|
---|
| 59 | #ifdef _DENSE_HASH_MAP_H_
|
---|
| 60 | /** Represents a google::dense_hash_map
|
---|
| 61 | *
|
---|
| 62 | * \see RandomSetter
|
---|
| 63 | */
|
---|
| 64 | template<typename Scalar> struct GoogleDenseHashMapTraits
|
---|
| 65 | {
|
---|
| 66 | typedef int KeyType;
|
---|
| 67 | typedef google::dense_hash_map<KeyType,Scalar> Type;
|
---|
| 68 | enum {
|
---|
| 69 | IsSorted = 0
|
---|
| 70 | };
|
---|
| 71 |
|
---|
| 72 | static void setInvalidKey(Type& map, const KeyType& k)
|
---|
| 73 | { map.set_empty_key(k); }
|
---|
| 74 | };
|
---|
| 75 | #endif
|
---|
| 76 |
|
---|
| 77 | #ifdef _SPARSE_HASH_MAP_H_
|
---|
| 78 | /** Represents a google::sparse_hash_map
|
---|
| 79 | *
|
---|
| 80 | * \see RandomSetter
|
---|
| 81 | */
|
---|
| 82 | template<typename Scalar> struct GoogleSparseHashMapTraits
|
---|
| 83 | {
|
---|
| 84 | typedef int KeyType;
|
---|
| 85 | typedef google::sparse_hash_map<KeyType,Scalar> Type;
|
---|
| 86 | enum {
|
---|
| 87 | IsSorted = 0
|
---|
| 88 | };
|
---|
| 89 |
|
---|
| 90 | static void setInvalidKey(Type&, const KeyType&) {}
|
---|
| 91 | };
|
---|
| 92 | #endif
|
---|
| 93 |
|
---|
| 94 | /** \class RandomSetter
|
---|
| 95 | *
|
---|
| 96 | * \brief The RandomSetter is a wrapper object allowing to set/update a sparse matrix with random access
|
---|
| 97 | *
|
---|
| 98 | * \param SparseMatrixType the type of the sparse matrix we are updating
|
---|
| 99 | * \param MapTraits a traits class representing the map implementation used for the temporary sparse storage.
|
---|
| 100 | * Its default value depends on the system.
|
---|
| 101 | * \param OuterPacketBits defines the number of rows (or columns) manage by a single map object
|
---|
| 102 | * as a power of two exponent.
|
---|
| 103 | *
|
---|
| 104 | * This class temporarily represents a sparse matrix object using a generic map implementation allowing for
|
---|
| 105 | * efficient random access. The conversion from the compressed representation to a hash_map object is performed
|
---|
| 106 | * in the RandomSetter constructor, while the sparse matrix is updated back at destruction time. This strategy
|
---|
| 107 | * suggest the use of nested blocks as in this example:
|
---|
| 108 | *
|
---|
| 109 | * \code
|
---|
| 110 | * SparseMatrix<double> m(rows,cols);
|
---|
| 111 | * {
|
---|
| 112 | * RandomSetter<SparseMatrix<double> > w(m);
|
---|
| 113 | * // don't use m but w instead with read/write random access to the coefficients:
|
---|
| 114 | * for(;;)
|
---|
| 115 | * w(rand(),rand()) = rand;
|
---|
| 116 | * }
|
---|
| 117 | * // when w is deleted, the data are copied back to m
|
---|
| 118 | * // and m is ready to use.
|
---|
| 119 | * \endcode
|
---|
| 120 | *
|
---|
| 121 | * Since hash_map objects are not fully sorted, representing a full matrix as a single hash_map would
|
---|
| 122 | * involve a big and costly sort to update the compressed matrix back. To overcome this issue, a RandomSetter
|
---|
| 123 | * use multiple hash_map, each representing 2^OuterPacketBits columns or rows according to the storage order.
|
---|
| 124 | * To reach optimal performance, this value should be adjusted according to the average number of nonzeros
|
---|
| 125 | * per rows/columns.
|
---|
| 126 | *
|
---|
| 127 | * The possible values for the template parameter MapTraits are:
|
---|
| 128 | * - \b StdMapTraits: corresponds to std::map. (does not perform very well)
|
---|
| 129 | * - \b GnuHashMapTraits: corresponds to __gnu_cxx::hash_map (available only with GCC)
|
---|
| 130 | * - \b GoogleDenseHashMapTraits: corresponds to google::dense_hash_map (best efficiency, reasonable memory consumption)
|
---|
| 131 | * - \b GoogleSparseHashMapTraits: corresponds to google::sparse_hash_map (best memory consumption, relatively good performance)
|
---|
| 132 | *
|
---|
| 133 | * The default map implementation depends on the availability, and the preferred order is:
|
---|
| 134 | * GoogleSparseHashMapTraits, GnuHashMapTraits, and finally StdMapTraits.
|
---|
| 135 | *
|
---|
| 136 | * For performance and memory consumption reasons it is highly recommended to use one of
|
---|
| 137 | * the Google's hash_map implementation. To enable the support for them, you have two options:
|
---|
| 138 | * - \#include <google/dense_hash_map> yourself \b before Eigen/Sparse header
|
---|
| 139 | * - define EIGEN_GOOGLEHASH_SUPPORT
|
---|
| 140 | * In the later case the inclusion of <google/dense_hash_map> is made for you.
|
---|
| 141 | *
|
---|
| 142 | * \see http://code.google.com/p/google-sparsehash/
|
---|
| 143 | */
|
---|
| 144 | template<typename SparseMatrixType,
|
---|
| 145 | template <typename T> class MapTraits =
|
---|
| 146 | #if defined _DENSE_HASH_MAP_H_
|
---|
| 147 | GoogleDenseHashMapTraits
|
---|
| 148 | #elif defined _HASH_MAP
|
---|
| 149 | GnuHashMapTraits
|
---|
| 150 | #else
|
---|
| 151 | StdMapTraits
|
---|
| 152 | #endif
|
---|
| 153 | ,int OuterPacketBits = 6>
|
---|
| 154 | class RandomSetter
|
---|
| 155 | {
|
---|
| 156 | typedef typename SparseMatrixType::Scalar Scalar;
|
---|
| 157 | typedef typename SparseMatrixType::Index Index;
|
---|
| 158 |
|
---|
| 159 | struct ScalarWrapper
|
---|
| 160 | {
|
---|
| 161 | ScalarWrapper() : value(0) {}
|
---|
| 162 | Scalar value;
|
---|
| 163 | };
|
---|
| 164 | typedef typename MapTraits<ScalarWrapper>::KeyType KeyType;
|
---|
| 165 | typedef typename MapTraits<ScalarWrapper>::Type HashMapType;
|
---|
| 166 | static const int OuterPacketMask = (1 << OuterPacketBits) - 1;
|
---|
| 167 | enum {
|
---|
| 168 | SwapStorage = 1 - MapTraits<ScalarWrapper>::IsSorted,
|
---|
| 169 | TargetRowMajor = (SparseMatrixType::Flags & RowMajorBit) ? 1 : 0,
|
---|
| 170 | SetterRowMajor = SwapStorage ? 1-TargetRowMajor : TargetRowMajor
|
---|
| 171 | };
|
---|
| 172 |
|
---|
| 173 | public:
|
---|
| 174 |
|
---|
| 175 | /** Constructs a random setter object from the sparse matrix \a target
|
---|
| 176 | *
|
---|
| 177 | * Note that the initial value of \a target are imported. If you want to re-set
|
---|
| 178 | * a sparse matrix from scratch, then you must set it to zero first using the
|
---|
| 179 | * setZero() function.
|
---|
| 180 | */
|
---|
| 181 | inline RandomSetter(SparseMatrixType& target)
|
---|
| 182 | : mp_target(&target)
|
---|
| 183 | {
|
---|
| 184 | const Index outerSize = SwapStorage ? target.innerSize() : target.outerSize();
|
---|
| 185 | const Index innerSize = SwapStorage ? target.outerSize() : target.innerSize();
|
---|
| 186 | m_outerPackets = outerSize >> OuterPacketBits;
|
---|
| 187 | if (outerSize&OuterPacketMask)
|
---|
| 188 | m_outerPackets += 1;
|
---|
| 189 | m_hashmaps = new HashMapType[m_outerPackets];
|
---|
| 190 | // compute number of bits needed to store inner indices
|
---|
| 191 | Index aux = innerSize - 1;
|
---|
| 192 | m_keyBitsOffset = 0;
|
---|
| 193 | while (aux)
|
---|
| 194 | {
|
---|
| 195 | ++m_keyBitsOffset;
|
---|
| 196 | aux = aux >> 1;
|
---|
| 197 | }
|
---|
| 198 | KeyType ik = (1<<(OuterPacketBits+m_keyBitsOffset));
|
---|
| 199 | for (Index k=0; k<m_outerPackets; ++k)
|
---|
| 200 | MapTraits<ScalarWrapper>::setInvalidKey(m_hashmaps[k],ik);
|
---|
| 201 |
|
---|
| 202 | // insert current coeffs
|
---|
| 203 | for (Index j=0; j<mp_target->outerSize(); ++j)
|
---|
| 204 | for (typename SparseMatrixType::InnerIterator it(*mp_target,j); it; ++it)
|
---|
| 205 | (*this)(TargetRowMajor?j:it.index(), TargetRowMajor?it.index():j) = it.value();
|
---|
| 206 | }
|
---|
| 207 |
|
---|
| 208 | /** Destructor updating back the sparse matrix target */
|
---|
| 209 | ~RandomSetter()
|
---|
| 210 | {
|
---|
| 211 | KeyType keyBitsMask = (1<<m_keyBitsOffset)-1;
|
---|
| 212 | if (!SwapStorage) // also means the map is sorted
|
---|
| 213 | {
|
---|
| 214 | mp_target->setZero();
|
---|
| 215 | mp_target->makeCompressed();
|
---|
| 216 | mp_target->reserve(nonZeros());
|
---|
| 217 | Index prevOuter = -1;
|
---|
| 218 | for (Index k=0; k<m_outerPackets; ++k)
|
---|
| 219 | {
|
---|
| 220 | const Index outerOffset = (1<<OuterPacketBits) * k;
|
---|
| 221 | typename HashMapType::iterator end = m_hashmaps[k].end();
|
---|
| 222 | for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it!=end; ++it)
|
---|
| 223 | {
|
---|
| 224 | const Index outer = (it->first >> m_keyBitsOffset) + outerOffset;
|
---|
| 225 | const Index inner = it->first & keyBitsMask;
|
---|
| 226 | if (prevOuter!=outer)
|
---|
| 227 | {
|
---|
| 228 | for (Index j=prevOuter+1;j<=outer;++j)
|
---|
| 229 | mp_target->startVec(j);
|
---|
| 230 | prevOuter = outer;
|
---|
| 231 | }
|
---|
| 232 | mp_target->insertBackByOuterInner(outer, inner) = it->second.value;
|
---|
| 233 | }
|
---|
| 234 | }
|
---|
| 235 | mp_target->finalize();
|
---|
| 236 | }
|
---|
| 237 | else
|
---|
| 238 | {
|
---|
| 239 | VectorXi positions(mp_target->outerSize());
|
---|
| 240 | positions.setZero();
|
---|
| 241 | // pass 1
|
---|
| 242 | for (Index k=0; k<m_outerPackets; ++k)
|
---|
| 243 | {
|
---|
| 244 | typename HashMapType::iterator end = m_hashmaps[k].end();
|
---|
| 245 | for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it!=end; ++it)
|
---|
| 246 | {
|
---|
| 247 | const Index outer = it->first & keyBitsMask;
|
---|
| 248 | ++positions[outer];
|
---|
| 249 | }
|
---|
| 250 | }
|
---|
| 251 | // prefix sum
|
---|
| 252 | Index count = 0;
|
---|
| 253 | for (Index j=0; j<mp_target->outerSize(); ++j)
|
---|
| 254 | {
|
---|
| 255 | Index tmp = positions[j];
|
---|
| 256 | mp_target->outerIndexPtr()[j] = count;
|
---|
| 257 | positions[j] = count;
|
---|
| 258 | count += tmp;
|
---|
| 259 | }
|
---|
| 260 | mp_target->makeCompressed();
|
---|
| 261 | mp_target->outerIndexPtr()[mp_target->outerSize()] = count;
|
---|
| 262 | mp_target->resizeNonZeros(count);
|
---|
| 263 | // pass 2
|
---|
| 264 | for (Index k=0; k<m_outerPackets; ++k)
|
---|
| 265 | {
|
---|
| 266 | const Index outerOffset = (1<<OuterPacketBits) * k;
|
---|
| 267 | typename HashMapType::iterator end = m_hashmaps[k].end();
|
---|
| 268 | for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it!=end; ++it)
|
---|
| 269 | {
|
---|
| 270 | const Index inner = (it->first >> m_keyBitsOffset) + outerOffset;
|
---|
| 271 | const Index outer = it->first & keyBitsMask;
|
---|
| 272 | // sorted insertion
|
---|
| 273 | // Note that we have to deal with at most 2^OuterPacketBits unsorted coefficients,
|
---|
| 274 | // moreover those 2^OuterPacketBits coeffs are likely to be sparse, an so only a
|
---|
| 275 | // small fraction of them have to be sorted, whence the following simple procedure:
|
---|
| 276 | Index posStart = mp_target->outerIndexPtr()[outer];
|
---|
| 277 | Index i = (positions[outer]++) - 1;
|
---|
| 278 | while ( (i >= posStart) && (mp_target->innerIndexPtr()[i] > inner) )
|
---|
| 279 | {
|
---|
| 280 | mp_target->valuePtr()[i+1] = mp_target->valuePtr()[i];
|
---|
| 281 | mp_target->innerIndexPtr()[i+1] = mp_target->innerIndexPtr()[i];
|
---|
| 282 | --i;
|
---|
| 283 | }
|
---|
| 284 | mp_target->innerIndexPtr()[i+1] = inner;
|
---|
| 285 | mp_target->valuePtr()[i+1] = it->second.value;
|
---|
| 286 | }
|
---|
| 287 | }
|
---|
| 288 | }
|
---|
| 289 | delete[] m_hashmaps;
|
---|
| 290 | }
|
---|
| 291 |
|
---|
| 292 | /** \returns a reference to the coefficient at given coordinates \a row, \a col */
|
---|
| 293 | Scalar& operator() (Index row, Index col)
|
---|
| 294 | {
|
---|
| 295 | const Index outer = SetterRowMajor ? row : col;
|
---|
| 296 | const Index inner = SetterRowMajor ? col : row;
|
---|
| 297 | const Index outerMajor = outer >> OuterPacketBits; // index of the packet/map
|
---|
| 298 | const Index outerMinor = outer & OuterPacketMask; // index of the inner vector in the packet
|
---|
| 299 | const KeyType key = (KeyType(outerMinor)<<m_keyBitsOffset) | inner;
|
---|
| 300 | return m_hashmaps[outerMajor][key].value;
|
---|
| 301 | }
|
---|
| 302 |
|
---|
| 303 | /** \returns the number of non zero coefficients
|
---|
| 304 | *
|
---|
| 305 | * \note According to the underlying map/hash_map implementation,
|
---|
| 306 | * this function might be quite expensive.
|
---|
| 307 | */
|
---|
| 308 | Index nonZeros() const
|
---|
| 309 | {
|
---|
| 310 | Index nz = 0;
|
---|
| 311 | for (Index k=0; k<m_outerPackets; ++k)
|
---|
| 312 | nz += static_cast<Index>(m_hashmaps[k].size());
|
---|
| 313 | return nz;
|
---|
| 314 | }
|
---|
| 315 |
|
---|
| 316 |
|
---|
| 317 | protected:
|
---|
| 318 |
|
---|
| 319 | HashMapType* m_hashmaps;
|
---|
| 320 | SparseMatrixType* mp_target;
|
---|
| 321 | Index m_outerPackets;
|
---|
| 322 | unsigned char m_keyBitsOffset;
|
---|
| 323 | };
|
---|
| 324 |
|
---|
| 325 | } // end namespace Eigen
|
---|
| 326 |
|
---|
| 327 | #endif // EIGEN_RANDOMSETTER_H
|
---|