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
|
---|