RandomSetter.h
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 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_RANDOMSETTER_H
26 #define EIGEN_RANDOMSETTER_H
27 
28 namespace Eigen {
29 
34 template<typename Scalar> struct StdMapTraits
35 {
36  typedef int KeyType;
37  typedef std::map<KeyType,Scalar> Type;
38  enum {
39  IsSorted = 1
40  };
41 
42  static void setInvalidKey(Type&, const KeyType&) {}
43 };
44 
45 #ifdef EIGEN_UNORDERED_MAP_SUPPORT
46 
62 template<typename Scalar> struct StdUnorderedMapTraits
63 {
64  typedef int KeyType;
65  typedef std::unordered_map<KeyType,Scalar> Type;
66  enum {
67  IsSorted = 0
68  };
69 
70  static void setInvalidKey(Type&, const KeyType&) {}
71 };
72 #endif // EIGEN_UNORDERED_MAP_SUPPORT
73 
74 #ifdef _DENSE_HASH_MAP_H_
75 
79 template<typename Scalar> struct GoogleDenseHashMapTraits
80 {
81  typedef int KeyType;
82  typedef google::dense_hash_map<KeyType,Scalar> Type;
83  enum {
84  IsSorted = 0
85  };
86 
87  static void setInvalidKey(Type& map, const KeyType& k)
88  { map.set_empty_key(k); }
89 };
90 #endif
91 
92 #ifdef _SPARSE_HASH_MAP_H_
93 
97 template<typename Scalar> struct GoogleSparseHashMapTraits
98 {
99  typedef int KeyType;
100  typedef google::sparse_hash_map<KeyType,Scalar> Type;
101  enum {
102  IsSorted = 0
103  };
104 
105  static void setInvalidKey(Type&, const KeyType&) {}
106 };
107 #endif
108 
159 template<typename SparseMatrixType,
160  template <typename T> class MapTraits =
161 #if defined _DENSE_HASH_MAP_H_
162  GoogleDenseHashMapTraits
163 #elif defined _HASH_MAP
164  GnuHashMapTraits
165 #else
166  StdMapTraits
167 #endif
168  ,int OuterPacketBits = 6>
170 {
171  typedef typename SparseMatrixType::Scalar Scalar;
172  typedef typename SparseMatrixType::Index Index;
173 
174  struct ScalarWrapper
175  {
176  ScalarWrapper() : value(0) {}
177  Scalar value;
178  };
179  typedef typename MapTraits<ScalarWrapper>::KeyType KeyType;
180  typedef typename MapTraits<ScalarWrapper>::Type HashMapType;
181  static const int OuterPacketMask = (1 << OuterPacketBits) - 1;
182  enum {
183  SwapStorage = 1 - MapTraits<ScalarWrapper>::IsSorted,
184  TargetRowMajor = (SparseMatrixType::Flags & RowMajorBit) ? 1 : 0,
185  SetterRowMajor = SwapStorage ? 1-TargetRowMajor : TargetRowMajor
186  };
187 
188  public:
189 
196  inline RandomSetter(SparseMatrixType& target)
197  : mp_target(&target)
198  {
199  const Index outerSize = SwapStorage ? target.innerSize() : target.outerSize();
200  const Index innerSize = SwapStorage ? target.outerSize() : target.innerSize();
201  m_outerPackets = outerSize >> OuterPacketBits;
202  if (outerSize&OuterPacketMask)
203  m_outerPackets += 1;
204  m_hashmaps = new HashMapType[m_outerPackets];
205  // compute number of bits needed to store inner indices
206  Index aux = innerSize - 1;
207  m_keyBitsOffset = 0;
208  while (aux)
209  {
210  ++m_keyBitsOffset;
211  aux = aux >> 1;
212  }
213  KeyType ik = (1<<(OuterPacketBits+m_keyBitsOffset));
214  for (Index k=0; k<m_outerPackets; ++k)
215  MapTraits<ScalarWrapper>::setInvalidKey(m_hashmaps[k],ik);
216 
217  // insert current coeffs
218  for (Index j=0; j<mp_target->outerSize(); ++j)
219  for (typename SparseMatrixType::InnerIterator it(*mp_target,j); it; ++it)
220  (*this)(TargetRowMajor?j:it.index(), TargetRowMajor?it.index():j) = it.value();
221  }
222 
225  {
226  KeyType keyBitsMask = (1<<m_keyBitsOffset)-1;
227  if (!SwapStorage) // also means the map is sorted
228  {
229  mp_target->setZero();
230  mp_target->makeCompressed();
231  mp_target->reserve(nonZeros());
232  Index prevOuter = -1;
233  for (Index k=0; k<m_outerPackets; ++k)
234  {
235  const Index outerOffset = (1<<OuterPacketBits) * k;
236  typename HashMapType::iterator end = m_hashmaps[k].end();
237  for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it!=end; ++it)
238  {
239  const Index outer = (it->first >> m_keyBitsOffset) + outerOffset;
240  const Index inner = it->first & keyBitsMask;
241  if (prevOuter!=outer)
242  {
243  for (Index j=prevOuter+1;j<=outer;++j)
244  mp_target->startVec(j);
245  prevOuter = outer;
246  }
247  mp_target->insertBackByOuterInner(outer, inner) = it->second.value;
248  }
249  }
250  mp_target->finalize();
251  }
252  else
253  {
254  VectorXi positions(mp_target->outerSize());
255  positions.setZero();
256  // pass 1
257  for (Index k=0; k<m_outerPackets; ++k)
258  {
259  typename HashMapType::iterator end = m_hashmaps[k].end();
260  for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it!=end; ++it)
261  {
262  const Index outer = it->first & keyBitsMask;
263  ++positions[outer];
264  }
265  }
266  // prefix sum
267  Index count = 0;
268  for (Index j=0; j<mp_target->outerSize(); ++j)
269  {
270  Index tmp = positions[j];
271  mp_target->outerIndexPtr()[j] = count;
272  positions[j] = count;
273  count += tmp;
274  }
275  mp_target->makeCompressed();
276  mp_target->outerIndexPtr()[mp_target->outerSize()] = count;
277  mp_target->resizeNonZeros(count);
278  // pass 2
279  for (Index k=0; k<m_outerPackets; ++k)
280  {
281  const Index outerOffset = (1<<OuterPacketBits) * k;
282  typename HashMapType::iterator end = m_hashmaps[k].end();
283  for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it!=end; ++it)
284  {
285  const Index inner = (it->first >> m_keyBitsOffset) + outerOffset;
286  const Index outer = it->first & keyBitsMask;
287  // sorted insertion
288  // Note that we have to deal with at most 2^OuterPacketBits unsorted coefficients,
289  // moreover those 2^OuterPacketBits coeffs are likely to be sparse, an so only a
290  // small fraction of them have to be sorted, whence the following simple procedure:
291  Index posStart = mp_target->outerIndexPtr()[outer];
292  Index i = (positions[outer]++) - 1;
293  while ( (i >= posStart) && (mp_target->innerIndexPtr()[i] > inner) )
294  {
295  mp_target->valuePtr()[i+1] = mp_target->valuePtr()[i];
296  mp_target->innerIndexPtr()[i+1] = mp_target->innerIndexPtr()[i];
297  --i;
298  }
299  mp_target->innerIndexPtr()[i+1] = inner;
300  mp_target->valuePtr()[i+1] = it->second.value;
301  }
302  }
303  }
304  delete[] m_hashmaps;
305  }
306 
308  Scalar& operator() (Index row, Index col)
309  {
310  const Index outer = SetterRowMajor ? row : col;
311  const Index inner = SetterRowMajor ? col : row;
312  const Index outerMajor = outer >> OuterPacketBits; // index of the packet/map
313  const Index outerMinor = outer & OuterPacketMask; // index of the inner vector in the packet
314  const KeyType key = (KeyType(outerMinor)<<m_keyBitsOffset) | inner;
315  return m_hashmaps[outerMajor][key].value;
316  }
317 
323  Index nonZeros() const
324  {
325  Index nz = 0;
326  for (Index k=0; k<m_outerPackets; ++k)
327  nz += static_cast<Index>(m_hashmaps[k].size());
328  return nz;
329  }
330 
331 
332  protected:
333 
334  HashMapType* m_hashmaps;
335  SparseMatrixType* mp_target;
336  Index m_outerPackets;
337  unsigned char m_keyBitsOffset;
338 };
339 
340 } // end namespace Eigen
341 
342 #endif // EIGEN_RANDOMSETTER_H