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 setInvalidKeyStdMapTraits27 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 setInvalidKeyStdUnorderedMapTraits55 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 setInvalidKeyGoogleDenseHashMapTraits72 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 setInvalidKeyGoogleSparseHashMapTraits90 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 * \tparam SparseMatrixType the type of the sparse matrix we are updating 99 * \tparam MapTraits a traits class representing the map implementation used for the temporary sparse storage. 100 * Its default value depends on the system. 101 * \tparam 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::StorageIndex StorageIndex; 158 159 struct ScalarWrapper 160 { ScalarWrapperScalarWrapper161 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 */ RandomSetter(SparseMatrixType & target)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 */ ~RandomSetter()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 */ operator()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 = internal::convert_index<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 */ nonZeros()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 328