Home | History | Annotate | Download | only in SparseExtra
      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 (at) 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
    328