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