RandomSetter.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_RANDOMSETTER_H
00026 #define EIGEN_RANDOMSETTER_H
00027 
00032 template<typename Scalar> struct StdMapTraits
00033 {
00034   typedef int KeyType;
00035   typedef std::map<KeyType,Scalar> Type;
00036   enum {
00037     IsSorted = 1
00038   };
00039 
00040   static void setInvalidKey(Type&, const KeyType&) {}
00041 };
00042 
00043 #ifdef EIGEN_UNORDERED_MAP_SUPPORT
00044 
00060 template<typename Scalar> struct StdUnorderedMapTraits
00061 {
00062   typedef int KeyType;
00063   typedef std::unordered_map<KeyType,Scalar> Type;
00064   enum {
00065     IsSorted = 0
00066   };
00067 
00068   static void setInvalidKey(Type&, const KeyType&) {}
00069 };
00070 #endif // EIGEN_UNORDERED_MAP_SUPPORT
00071 
00072 #ifdef _DENSE_HASH_MAP_H_
00073 
00077 template<typename Scalar> struct GoogleDenseHashMapTraits
00078 {
00079   typedef int KeyType;
00080   typedef google::dense_hash_map<KeyType,Scalar> Type;
00081   enum {
00082     IsSorted = 0
00083   };
00084 
00085   static void setInvalidKey(Type& map, const KeyType& k)
00086   { map.set_empty_key(k); }
00087 };
00088 #endif
00089 
00090 #ifdef _SPARSE_HASH_MAP_H_
00091 
00095 template<typename Scalar> struct GoogleSparseHashMapTraits
00096 {
00097   typedef int KeyType;
00098   typedef google::sparse_hash_map<KeyType,Scalar> Type;
00099   enum {
00100     IsSorted = 0
00101   };
00102 
00103   static void setInvalidKey(Type&, const KeyType&) {}
00104 };
00105 #endif
00106 
00157 template<typename SparseMatrixType,
00158          template <typename T> class MapTraits =
00159 #if defined _DENSE_HASH_MAP_H_
00160           GoogleDenseHashMapTraits
00161 #elif defined _HASH_MAP
00162           GnuHashMapTraits
00163 #else
00164           StdMapTraits
00165 #endif
00166          ,int OuterPacketBits = 6>
00167 class RandomSetter
00168 {
00169     typedef typename SparseMatrixType::Scalar Scalar;
00170     typedef typename SparseMatrixType::Index Index;
00171 
00172     struct ScalarWrapper
00173     {
00174       ScalarWrapper() : value(0) {}
00175       Scalar value;
00176     };
00177     typedef typename MapTraits<ScalarWrapper>::KeyType KeyType;
00178     typedef typename MapTraits<ScalarWrapper>::Type HashMapType;
00179     static const int OuterPacketMask = (1 << OuterPacketBits) - 1;
00180     enum {
00181       SwapStorage = 1 - MapTraits<ScalarWrapper>::IsSorted,
00182       TargetRowMajor = (SparseMatrixType::Flags & RowMajorBit) ? 1 : 0,
00183       SetterRowMajor = SwapStorage ? 1-TargetRowMajor : TargetRowMajor,
00184       IsUpper = SparseMatrixType::Flags & Upper,
00185       IsLower = SparseMatrixType::Flags & Lower
00186     };
00187 
00188   public:
00189 
00196     inline RandomSetter(SparseMatrixType& target)
00197       : mp_target(&target)
00198     {
00199       const Index outerSize = SwapStorage ? target.innerSize() : target.outerSize();
00200       const Index innerSize = SwapStorage ? target.outerSize() : target.innerSize();
00201       m_outerPackets = outerSize >> OuterPacketBits;
00202       if (outerSize&OuterPacketMask)
00203         m_outerPackets += 1;
00204       m_hashmaps = new HashMapType[m_outerPackets];
00205       // compute number of bits needed to store inner indices
00206       Index aux = innerSize - 1;
00207       m_keyBitsOffset = 0;
00208       while (aux)
00209       {
00210         ++m_keyBitsOffset;
00211         aux = aux >> 1;
00212       }
00213       KeyType ik = (1<<(OuterPacketBits+m_keyBitsOffset));
00214       for (Index k=0; k<m_outerPackets; ++k)
00215         MapTraits<ScalarWrapper>::setInvalidKey(m_hashmaps[k],ik);
00216 
00217       // insert current coeffs
00218       for (Index j=0; j<mp_target->outerSize(); ++j)
00219         for (typename SparseMatrixType::InnerIterator it(*mp_target,j); it; ++it)
00220           (*this)(TargetRowMajor?j:it.index(), TargetRowMajor?it.index():j) = it.value();
00221     }
00222 
00224     ~RandomSetter()
00225     {
00226       KeyType keyBitsMask = (1<<m_keyBitsOffset)-1;
00227       if (!SwapStorage) // also means the map is sorted
00228       {
00229         mp_target->setZero();
00230         mp_target->reserve(nonZeros());
00231         Index prevOuter = -1;
00232         for (Index k=0; k<m_outerPackets; ++k)
00233         {
00234           const Index outerOffset = (1<<OuterPacketBits) * k;
00235           typename HashMapType::iterator end = m_hashmaps[k].end();
00236           for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it!=end; ++it)
00237           {
00238             const Index outer = (it->first >> m_keyBitsOffset) + outerOffset;
00239             const Index inner = it->first & keyBitsMask;
00240             if (prevOuter!=outer)
00241             {
00242               for (Index j=prevOuter+1;j<=outer;++j)
00243                 mp_target->startVec(j);
00244               prevOuter = outer;
00245             }
00246             mp_target->insertBackByOuterInner(outer, inner) = it->second.value;
00247           }
00248         }
00249         mp_target->finalize();
00250       }
00251       else
00252       {
00253         VectorXi positions(mp_target->outerSize());
00254         positions.setZero();
00255         // pass 1
00256         for (Index k=0; k<m_outerPackets; ++k)
00257         {
00258           typename HashMapType::iterator end = m_hashmaps[k].end();
00259           for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it!=end; ++it)
00260           {
00261             const Index outer = it->first & keyBitsMask;
00262             ++positions[outer];
00263           }
00264         }
00265         // prefix sum
00266         Index count = 0;
00267         for (Index j=0; j<mp_target->outerSize(); ++j)
00268         {
00269           Index tmp = positions[j];
00270           mp_target->_outerIndexPtr()[j] = count;
00271           positions[j] = count;
00272           count += tmp;
00273         }
00274         mp_target->_outerIndexPtr()[mp_target->outerSize()] = count;
00275         mp_target->resizeNonZeros(count);
00276         // pass 2
00277         for (Index k=0; k<m_outerPackets; ++k)
00278         {
00279           const Index outerOffset = (1<<OuterPacketBits) * k;
00280           typename HashMapType::iterator end = m_hashmaps[k].end();
00281           for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it!=end; ++it)
00282           {
00283             const Index inner = (it->first >> m_keyBitsOffset) + outerOffset;
00284             const Index outer = it->first & keyBitsMask;
00285             // sorted insertion
00286             // Note that we have to deal with at most 2^OuterPacketBits unsorted coefficients,
00287             // moreover those 2^OuterPacketBits coeffs are likely to be sparse, an so only a
00288             // small fraction of them have to be sorted, whence the following simple procedure:
00289             Index posStart = mp_target->_outerIndexPtr()[outer];
00290             Index i = (positions[outer]++) - 1;
00291             while ( (i >= posStart) && (mp_target->_innerIndexPtr()[i] > inner) )
00292             {
00293               mp_target->_valuePtr()[i+1] = mp_target->_valuePtr()[i];
00294               mp_target->_innerIndexPtr()[i+1] = mp_target->_innerIndexPtr()[i];
00295               --i;
00296             }
00297             mp_target->_innerIndexPtr()[i+1] = inner;
00298             mp_target->_valuePtr()[i+1] = it->second.value;
00299           }
00300         }
00301       }
00302       delete[] m_hashmaps;
00303     }
00304 
00306     Scalar& operator() (Index row, Index col)
00307     {
00308       eigen_assert(((!IsUpper) || (row<=col)) && "Invalid access to an upper triangular matrix");
00309       eigen_assert(((!IsLower) || (col<=row)) && "Invalid access to an upper triangular matrix");
00310       const Index outer = SetterRowMajor ? row : col;
00311       const Index inner = SetterRowMajor ? col : row;
00312       const Index outerMajor = outer >> OuterPacketBits; // index of the packet/map
00313       const Index outerMinor = outer & OuterPacketMask;  // index of the inner vector in the packet
00314       const KeyType key = (KeyType(outerMinor)<<m_keyBitsOffset) | inner;
00315       return m_hashmaps[outerMajor][key].value;
00316     }
00317 
00323     Index nonZeros() const
00324     {
00325       Index nz = 0;
00326       for (Index k=0; k<m_outerPackets; ++k)
00327         nz += static_cast<Index>(m_hashmaps[k].size());
00328       return nz;
00329     }
00330 
00331 
00332   protected:
00333 
00334     HashMapType* m_hashmaps;
00335     SparseMatrixType* mp_target;
00336     Index m_outerPackets;
00337     unsigned char m_keyBitsOffset;
00338 };
00339 
00340 #endif // EIGEN_RANDOMSETTER_H


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:16