$treeview $search $mathjax
Eigen-unsupported
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_RANDOMSETTER_H 00011 #define EIGEN_RANDOMSETTER_H 00012 00013 namespace Eigen { 00014 00019 template<typename Scalar> struct StdMapTraits 00020 { 00021 typedef int KeyType; 00022 typedef std::map<KeyType,Scalar> Type; 00023 enum { 00024 IsSorted = 1 00025 }; 00026 00027 static void setInvalidKey(Type&, const KeyType&) {} 00028 }; 00029 00030 #ifdef EIGEN_UNORDERED_MAP_SUPPORT 00031 00047 template<typename Scalar> struct StdUnorderedMapTraits 00048 { 00049 typedef int KeyType; 00050 typedef std::unordered_map<KeyType,Scalar> Type; 00051 enum { 00052 IsSorted = 0 00053 }; 00054 00055 static void setInvalidKey(Type&, const KeyType&) {} 00056 }; 00057 #endif // EIGEN_UNORDERED_MAP_SUPPORT 00058 00059 #ifdef _DENSE_HASH_MAP_H_ 00060 00064 template<typename Scalar> struct GoogleDenseHashMapTraits 00065 { 00066 typedef int KeyType; 00067 typedef google::dense_hash_map<KeyType,Scalar> Type; 00068 enum { 00069 IsSorted = 0 00070 }; 00071 00072 static void setInvalidKey(Type& map, const KeyType& k) 00073 { map.set_empty_key(k); } 00074 }; 00075 #endif 00076 00077 #ifdef _SPARSE_HASH_MAP_H_ 00078 00082 template<typename Scalar> struct GoogleSparseHashMapTraits 00083 { 00084 typedef int KeyType; 00085 typedef google::sparse_hash_map<KeyType,Scalar> Type; 00086 enum { 00087 IsSorted = 0 00088 }; 00089 00090 static void setInvalidKey(Type&, const KeyType&) {} 00091 }; 00092 #endif 00093 00144 template<typename SparseMatrixType, 00145 template <typename T> class MapTraits = 00146 #if defined _DENSE_HASH_MAP_H_ 00147 GoogleDenseHashMapTraits 00148 #elif defined _HASH_MAP 00149 GnuHashMapTraits 00150 #else 00151 StdMapTraits 00152 #endif 00153 ,int OuterPacketBits = 6> 00154 class RandomSetter 00155 { 00156 typedef typename SparseMatrixType::Scalar Scalar; 00157 typedef typename SparseMatrixType::Index Index; 00158 00159 struct ScalarWrapper 00160 { 00161 ScalarWrapper() : value(0) {} 00162 Scalar value; 00163 }; 00164 typedef typename MapTraits<ScalarWrapper>::KeyType KeyType; 00165 typedef typename MapTraits<ScalarWrapper>::Type HashMapType; 00166 static const int OuterPacketMask = (1 << OuterPacketBits) - 1; 00167 enum { 00168 SwapStorage = 1 - MapTraits<ScalarWrapper>::IsSorted, 00169 TargetRowMajor = (SparseMatrixType::Flags & RowMajorBit) ? 1 : 0, 00170 SetterRowMajor = SwapStorage ? 1-TargetRowMajor : TargetRowMajor 00171 }; 00172 00173 public: 00174 00181 inline RandomSetter(SparseMatrixType& target) 00182 : mp_target(&target) 00183 { 00184 const Index outerSize = SwapStorage ? target.innerSize() : target.outerSize(); 00185 const Index innerSize = SwapStorage ? target.outerSize() : target.innerSize(); 00186 m_outerPackets = outerSize >> OuterPacketBits; 00187 if (outerSize&OuterPacketMask) 00188 m_outerPackets += 1; 00189 m_hashmaps = new HashMapType[m_outerPackets]; 00190 // compute number of bits needed to store inner indices 00191 Index aux = innerSize - 1; 00192 m_keyBitsOffset = 0; 00193 while (aux) 00194 { 00195 ++m_keyBitsOffset; 00196 aux = aux >> 1; 00197 } 00198 KeyType ik = (1<<(OuterPacketBits+m_keyBitsOffset)); 00199 for (Index k=0; k<m_outerPackets; ++k) 00200 MapTraits<ScalarWrapper>::setInvalidKey(m_hashmaps[k],ik); 00201 00202 // insert current coeffs 00203 for (Index j=0; j<mp_target->outerSize(); ++j) 00204 for (typename SparseMatrixType::InnerIterator it(*mp_target,j); it; ++it) 00205 (*this)(TargetRowMajor?j:it.index(), TargetRowMajor?it.index():j) = it.value(); 00206 } 00207 00209 ~RandomSetter() 00210 { 00211 KeyType keyBitsMask = (1<<m_keyBitsOffset)-1; 00212 if (!SwapStorage) // also means the map is sorted 00213 { 00214 mp_target->setZero(); 00215 mp_target->makeCompressed(); 00216 mp_target->reserve(nonZeros()); 00217 Index prevOuter = -1; 00218 for (Index k=0; k<m_outerPackets; ++k) 00219 { 00220 const Index outerOffset = (1<<OuterPacketBits) * k; 00221 typename HashMapType::iterator end = m_hashmaps[k].end(); 00222 for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it!=end; ++it) 00223 { 00224 const Index outer = (it->first >> m_keyBitsOffset) + outerOffset; 00225 const Index inner = it->first & keyBitsMask; 00226 if (prevOuter!=outer) 00227 { 00228 for (Index j=prevOuter+1;j<=outer;++j) 00229 mp_target->startVec(j); 00230 prevOuter = outer; 00231 } 00232 mp_target->insertBackByOuterInner(outer, inner) = it->second.value; 00233 } 00234 } 00235 mp_target->finalize(); 00236 } 00237 else 00238 { 00239 VectorXi positions(mp_target->outerSize()); 00240 positions.setZero(); 00241 // pass 1 00242 for (Index k=0; k<m_outerPackets; ++k) 00243 { 00244 typename HashMapType::iterator end = m_hashmaps[k].end(); 00245 for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it!=end; ++it) 00246 { 00247 const Index outer = it->first & keyBitsMask; 00248 ++positions[outer]; 00249 } 00250 } 00251 // prefix sum 00252 Index count = 0; 00253 for (Index j=0; j<mp_target->outerSize(); ++j) 00254 { 00255 Index tmp = positions[j]; 00256 mp_target->outerIndexPtr()[j] = count; 00257 positions[j] = count; 00258 count += tmp; 00259 } 00260 mp_target->makeCompressed(); 00261 mp_target->outerIndexPtr()[mp_target->outerSize()] = count; 00262 mp_target->resizeNonZeros(count); 00263 // pass 2 00264 for (Index k=0; k<m_outerPackets; ++k) 00265 { 00266 const Index outerOffset = (1<<OuterPacketBits) * k; 00267 typename HashMapType::iterator end = m_hashmaps[k].end(); 00268 for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it!=end; ++it) 00269 { 00270 const Index inner = (it->first >> m_keyBitsOffset) + outerOffset; 00271 const Index outer = it->first & keyBitsMask; 00272 // sorted insertion 00273 // Note that we have to deal with at most 2^OuterPacketBits unsorted coefficients, 00274 // moreover those 2^OuterPacketBits coeffs are likely to be sparse, an so only a 00275 // small fraction of them have to be sorted, whence the following simple procedure: 00276 Index posStart = mp_target->outerIndexPtr()[outer]; 00277 Index i = (positions[outer]++) - 1; 00278 while ( (i >= posStart) && (mp_target->innerIndexPtr()[i] > inner) ) 00279 { 00280 mp_target->valuePtr()[i+1] = mp_target->valuePtr()[i]; 00281 mp_target->innerIndexPtr()[i+1] = mp_target->innerIndexPtr()[i]; 00282 --i; 00283 } 00284 mp_target->innerIndexPtr()[i+1] = inner; 00285 mp_target->valuePtr()[i+1] = it->second.value; 00286 } 00287 } 00288 } 00289 delete[] m_hashmaps; 00290 } 00291 00293 Scalar& operator() (Index row, Index col) 00294 { 00295 const Index outer = SetterRowMajor ? row : col; 00296 const Index inner = SetterRowMajor ? col : row; 00297 const Index outerMajor = outer >> OuterPacketBits; // index of the packet/map 00298 const Index outerMinor = outer & OuterPacketMask; // index of the inner vector in the packet 00299 const KeyType key = (KeyType(outerMinor)<<m_keyBitsOffset) | inner; 00300 return m_hashmaps[outerMajor][key].value; 00301 } 00302 00308 Index nonZeros() const 00309 { 00310 Index nz = 0; 00311 for (Index k=0; k<m_outerPackets; ++k) 00312 nz += static_cast<Index>(m_hashmaps[k].size()); 00313 return nz; 00314 } 00315 00316 00317 protected: 00318 00319 HashMapType* m_hashmaps; 00320 SparseMatrixType* mp_target; 00321 Index m_outerPackets; 00322 unsigned char m_keyBitsOffset; 00323 }; 00324 00325 } // end namespace Eigen 00326 00327 #endif // EIGEN_RANDOMSETTER_H