Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:57:08

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009 
0010 #ifndef EIGEN_RANDOMSETTER_H
0011 #define EIGEN_RANDOMSETTER_H
0012 
0013 #if defined(EIGEN_GOOGLEHASH_SUPPORT)
0014 // Ensure the ::google namespace exists, required for checking existence of 
0015 // ::google::dense_hash_map and ::google::sparse_hash_map.
0016 namespace google {}
0017 #endif
0018 
0019 namespace Eigen {
0020 
0021 /** Represents a std::map
0022   *
0023   * \see RandomSetter
0024   */
0025 template<typename Scalar> struct StdMapTraits
0026 {
0027   typedef int KeyType;
0028   typedef std::map<KeyType,Scalar> Type;
0029   enum {
0030     IsSorted = 1
0031   };
0032 
0033   static void setInvalidKey(Type&, const KeyType&) {}
0034 };
0035 
0036 #ifdef EIGEN_UNORDERED_MAP_SUPPORT
0037 /** Represents a std::unordered_map
0038   *
0039   * To use it you need to both define EIGEN_UNORDERED_MAP_SUPPORT and include the unordered_map header file
0040   * yourself making sure that unordered_map is defined in the std namespace.
0041   *
0042   * For instance, with current version of gcc you can either enable C++0x standard (-std=c++0x) or do:
0043   * \code
0044   * #include <tr1/unordered_map>
0045   * #define EIGEN_UNORDERED_MAP_SUPPORT
0046   * namespace std {
0047   *   using std::tr1::unordered_map;
0048   * }
0049   * \endcode
0050   *
0051   * \see RandomSetter
0052   */
0053 template<typename Scalar> struct StdUnorderedMapTraits
0054 {
0055   typedef int KeyType;
0056   typedef std::unordered_map<KeyType,Scalar> Type;
0057   enum {
0058     IsSorted = 0
0059   };
0060 
0061   static void setInvalidKey(Type&, const KeyType&) {}
0062 };
0063 #endif // EIGEN_UNORDERED_MAP_SUPPORT
0064 
0065 #if defined(EIGEN_GOOGLEHASH_SUPPORT)
0066 
0067 namespace google {
0068   
0069 // Namespace work-around, since sometimes dense_hash_map and sparse_hash_map
0070 // are in the global namespace, and other times they are under ::google.
0071 using namespace ::google;
0072 
0073 template<typename KeyType, typename Scalar>
0074 struct DenseHashMap {
0075   typedef dense_hash_map<KeyType, Scalar> type;
0076 };
0077 
0078 template<typename KeyType, typename Scalar>
0079 struct SparseHashMap {
0080   typedef sparse_hash_map<KeyType, Scalar> type;
0081 };
0082 
0083 } // namespace google
0084 
0085 /** Represents a google::dense_hash_map
0086   *
0087   * \see RandomSetter
0088   */
0089 template<typename Scalar> struct GoogleDenseHashMapTraits
0090 {
0091   typedef int KeyType;
0092   typedef typename google::DenseHashMap<KeyType,Scalar>::type Type;
0093   enum {
0094     IsSorted = 0
0095   };
0096 
0097   static void setInvalidKey(Type& map, const KeyType& k)
0098   { map.set_empty_key(k); }
0099 };
0100 
0101 /** Represents a google::sparse_hash_map
0102   *
0103   * \see RandomSetter
0104   */
0105 template<typename Scalar> struct GoogleSparseHashMapTraits
0106 {
0107   typedef int KeyType;
0108   typedef typename google::SparseHashMap<KeyType,Scalar>::type Type;
0109   enum {
0110     IsSorted = 0
0111   };
0112 
0113   static void setInvalidKey(Type&, const KeyType&) {}
0114 };
0115 #endif
0116 
0117 /** \class RandomSetter
0118   *
0119   * \brief The RandomSetter is a wrapper object allowing to set/update a sparse matrix with random access
0120   *
0121   * \tparam SparseMatrixType the type of the sparse matrix we are updating
0122   * \tparam MapTraits a traits class representing the map implementation used for the temporary sparse storage.
0123   *                  Its default value depends on the system.
0124   * \tparam OuterPacketBits defines the number of rows (or columns) manage by a single map object
0125   *                        as a power of two exponent.
0126   *
0127   * This class temporarily represents a sparse matrix object using a generic map implementation allowing for
0128   * efficient random access. The conversion from the compressed representation to a hash_map object is performed
0129   * in the RandomSetter constructor, while the sparse matrix is updated back at destruction time. This strategy
0130   * suggest the use of nested blocks as in this example:
0131   *
0132   * \code
0133   * SparseMatrix<double> m(rows,cols);
0134   * {
0135   *   RandomSetter<SparseMatrix<double> > w(m);
0136   *   // don't use m but w instead with read/write random access to the coefficients:
0137   *   for(;;)
0138   *     w(rand(),rand()) = rand;
0139   * }
0140   * // when w is deleted, the data are copied back to m
0141   * // and m is ready to use.
0142   * \endcode
0143   *
0144   * Since hash_map objects are not fully sorted, representing a full matrix as a single hash_map would
0145   * involve a big and costly sort to update the compressed matrix back. To overcome this issue, a RandomSetter
0146   * use multiple hash_map, each representing 2^OuterPacketBits columns or rows according to the storage order.
0147   * To reach optimal performance, this value should be adjusted according to the average number of nonzeros
0148   * per rows/columns.
0149   *
0150   * The possible values for the template parameter MapTraits are:
0151   *  - \b StdMapTraits: corresponds to std::map. (does not perform very well)
0152   *  - \b GnuHashMapTraits: corresponds to __gnu_cxx::hash_map (available only with GCC)
0153   *  - \b GoogleDenseHashMapTraits: corresponds to google::dense_hash_map (best efficiency, reasonable memory consumption)
0154   *  - \b GoogleSparseHashMapTraits: corresponds to google::sparse_hash_map (best memory consumption, relatively good performance)
0155   *
0156   * The default map implementation depends on the availability, and the preferred order is:
0157   * GoogleSparseHashMapTraits, GnuHashMapTraits, and finally StdMapTraits.
0158   *
0159   * For performance and memory consumption reasons it is highly recommended to use one of
0160   * Google's hash_map implementations. To enable the support for them, you must define
0161   * EIGEN_GOOGLEHASH_SUPPORT. This will include both <google/dense_hash_map> and
0162   * <google/sparse_hash_map> for you.
0163   *
0164   * \see https://github.com/sparsehash/sparsehash
0165   */
0166 template<typename SparseMatrixType,
0167          template <typename T> class MapTraits =
0168 #if defined(EIGEN_GOOGLEHASH_SUPPORT)
0169           GoogleDenseHashMapTraits
0170 #elif defined(_HASH_MAP)
0171           GnuHashMapTraits
0172 #else
0173           StdMapTraits
0174 #endif
0175          ,int OuterPacketBits = 6>
0176 class RandomSetter
0177 {
0178     typedef typename SparseMatrixType::Scalar Scalar;
0179     typedef typename SparseMatrixType::StorageIndex StorageIndex;
0180 
0181     struct ScalarWrapper
0182     {
0183       ScalarWrapper() : value(0) {}
0184       Scalar value;
0185     };
0186     typedef typename MapTraits<ScalarWrapper>::KeyType KeyType;
0187     typedef typename MapTraits<ScalarWrapper>::Type HashMapType;
0188     static const int OuterPacketMask = (1 << OuterPacketBits) - 1;
0189     enum {
0190       SwapStorage = 1 - MapTraits<ScalarWrapper>::IsSorted,
0191       TargetRowMajor = (SparseMatrixType::Flags & RowMajorBit) ? 1 : 0,
0192       SetterRowMajor = SwapStorage ? 1-TargetRowMajor : TargetRowMajor
0193     };
0194 
0195   public:
0196 
0197     /** Constructs a random setter object from the sparse matrix \a target
0198       *
0199       * Note that the initial value of \a target are imported. If you want to re-set
0200       * a sparse matrix from scratch, then you must set it to zero first using the
0201       * setZero() function.
0202       */
0203     inline RandomSetter(SparseMatrixType& target)
0204       : mp_target(&target)
0205     {
0206       const Index outerSize = SwapStorage ? target.innerSize() : target.outerSize();
0207       const Index innerSize = SwapStorage ? target.outerSize() : target.innerSize();
0208       m_outerPackets = outerSize >> OuterPacketBits;
0209       if (outerSize&OuterPacketMask)
0210         m_outerPackets += 1;
0211       m_hashmaps = new HashMapType[m_outerPackets];
0212       // compute number of bits needed to store inner indices
0213       Index aux = innerSize - 1;
0214       m_keyBitsOffset = 0;
0215       while (aux)
0216       {
0217         ++m_keyBitsOffset;
0218         aux = aux >> 1;
0219       }
0220       KeyType ik = (1<<(OuterPacketBits+m_keyBitsOffset));
0221       for (Index k=0; k<m_outerPackets; ++k)
0222         MapTraits<ScalarWrapper>::setInvalidKey(m_hashmaps[k],ik);
0223 
0224       // insert current coeffs
0225       for (Index j=0; j<mp_target->outerSize(); ++j)
0226         for (typename SparseMatrixType::InnerIterator it(*mp_target,j); it; ++it)
0227           (*this)(TargetRowMajor?j:it.index(), TargetRowMajor?it.index():j) = it.value();
0228     }
0229 
0230     /** Destructor updating back the sparse matrix target */
0231     ~RandomSetter()
0232     {
0233       KeyType keyBitsMask = (1<<m_keyBitsOffset)-1;
0234       if (!SwapStorage) // also means the map is sorted
0235       {
0236         mp_target->setZero();
0237         mp_target->makeCompressed();
0238         mp_target->reserve(nonZeros());
0239         Index prevOuter = -1;
0240         for (Index k=0; k<m_outerPackets; ++k)
0241         {
0242           const Index outerOffset = (1<<OuterPacketBits) * k;
0243           typename HashMapType::iterator end = m_hashmaps[k].end();
0244           for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it!=end; ++it)
0245           {
0246             const Index outer = (it->first >> m_keyBitsOffset) + outerOffset;
0247             const Index inner = it->first & keyBitsMask;
0248             if (prevOuter!=outer)
0249             {
0250               for (Index j=prevOuter+1;j<=outer;++j)
0251                 mp_target->startVec(j);
0252               prevOuter = outer;
0253             }
0254             mp_target->insertBackByOuterInner(outer, inner) = it->second.value;
0255           }
0256         }
0257         mp_target->finalize();
0258       }
0259       else
0260       {
0261         VectorXi positions(mp_target->outerSize());
0262         positions.setZero();
0263         // pass 1
0264         for (Index k=0; k<m_outerPackets; ++k)
0265         {
0266           typename HashMapType::iterator end = m_hashmaps[k].end();
0267           for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it!=end; ++it)
0268           {
0269             const Index outer = it->first & keyBitsMask;
0270             ++positions[outer];
0271           }
0272         }
0273         // prefix sum
0274         StorageIndex count = 0;
0275         for (Index j=0; j<mp_target->outerSize(); ++j)
0276         {
0277           StorageIndex tmp = positions[j];
0278           mp_target->outerIndexPtr()[j] = count;
0279           positions[j] = count;
0280           count += tmp;
0281         }
0282         mp_target->makeCompressed();
0283         mp_target->outerIndexPtr()[mp_target->outerSize()] = count;
0284         mp_target->resizeNonZeros(count);
0285         // pass 2
0286         for (Index k=0; k<m_outerPackets; ++k)
0287         {
0288           const Index outerOffset = (1<<OuterPacketBits) * k;
0289           typename HashMapType::iterator end = m_hashmaps[k].end();
0290           for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it!=end; ++it)
0291           {
0292             const Index inner = (it->first >> m_keyBitsOffset) + outerOffset;
0293             const Index outer = it->first & keyBitsMask;
0294             // sorted insertion
0295             // Note that we have to deal with at most 2^OuterPacketBits unsorted coefficients,
0296             // moreover those 2^OuterPacketBits coeffs are likely to be sparse, an so only a
0297             // small fraction of them have to be sorted, whence the following simple procedure:
0298             Index posStart = mp_target->outerIndexPtr()[outer];
0299             Index i = (positions[outer]++) - 1;
0300             while ( (i >= posStart) && (mp_target->innerIndexPtr()[i] > inner) )
0301             {
0302               mp_target->valuePtr()[i+1] = mp_target->valuePtr()[i];
0303               mp_target->innerIndexPtr()[i+1] = mp_target->innerIndexPtr()[i];
0304               --i;
0305             }
0306             mp_target->innerIndexPtr()[i+1] = internal::convert_index<StorageIndex>(inner);
0307             mp_target->valuePtr()[i+1] = it->second.value;
0308           }
0309         }
0310       }
0311       delete[] m_hashmaps;
0312     }
0313 
0314     /** \returns a reference to the coefficient at given coordinates \a row, \a col */
0315     Scalar& operator() (Index row, Index col)
0316     {
0317       const Index outer = SetterRowMajor ? row : col;
0318       const Index inner = SetterRowMajor ? col : row;
0319       const Index outerMajor = outer >> OuterPacketBits; // index of the packet/map
0320       const Index outerMinor = outer & OuterPacketMask;  // index of the inner vector in the packet
0321       const KeyType key = internal::convert_index<KeyType>((outerMinor<<m_keyBitsOffset) | inner);
0322       return m_hashmaps[outerMajor][key].value;
0323     }
0324 
0325     /** \returns the number of non zero coefficients
0326       *
0327       * \note According to the underlying map/hash_map implementation,
0328       * this function might be quite expensive.
0329       */
0330     Index nonZeros() const
0331     {
0332       Index nz = 0;
0333       for (Index k=0; k<m_outerPackets; ++k)
0334         nz += static_cast<Index>(m_hashmaps[k].size());
0335       return nz;
0336     }
0337 
0338 
0339   protected:
0340 
0341     HashMapType* m_hashmaps;
0342     SparseMatrixType* mp_target;
0343     Index m_outerPackets;
0344     unsigned char m_keyBitsOffset;
0345 };
0346 
0347 } // end namespace Eigen
0348 
0349 #endif // EIGEN_RANDOMSETTER_H