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@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_AMBIVECTOR_H 11#define EIGEN_AMBIVECTOR_H 12 13namespace Eigen { 14 15namespace internal { 16 17/** \internal 18 * Hybrid sparse/dense vector class designed for intensive read-write operations. 19 * 20 * See BasicSparseLLT and SparseProduct for usage examples. 21 */ 22template<typename _Scalar, typename _StorageIndex> 23class AmbiVector 24{ 25 public: 26 typedef _Scalar Scalar; 27 typedef _StorageIndex StorageIndex; 28 typedef typename NumTraits<Scalar>::Real RealScalar; 29 30 explicit AmbiVector(Index size) 31 : m_buffer(0), m_zero(0), m_size(0), m_allocatedSize(0), m_allocatedElements(0), m_mode(-1) 32 { 33 resize(size); 34 } 35 36 void init(double estimatedDensity); 37 void init(int mode); 38 39 Index nonZeros() const; 40 41 /** Specifies a sub-vector to work on */ 42 void setBounds(Index start, Index end) { m_start = convert_index(start); m_end = convert_index(end); } 43 44 void setZero(); 45 46 void restart(); 47 Scalar& coeffRef(Index i); 48 Scalar& coeff(Index i); 49 50 class Iterator; 51 52 ~AmbiVector() { delete[] m_buffer; } 53 54 void resize(Index size) 55 { 56 if (m_allocatedSize < size) 57 reallocate(size); 58 m_size = convert_index(size); 59 } 60 61 StorageIndex size() const { return m_size; } 62 63 protected: 64 StorageIndex convert_index(Index idx) 65 { 66 return internal::convert_index<StorageIndex>(idx); 67 } 68 69 void reallocate(Index size) 70 { 71 // if the size of the matrix is not too large, let's allocate a bit more than needed such 72 // that we can handle dense vector even in sparse mode. 73 delete[] m_buffer; 74 if (size<1000) 75 { 76 Index allocSize = (size * sizeof(ListEl) + sizeof(Scalar) - 1)/sizeof(Scalar); 77 m_allocatedElements = convert_index((allocSize*sizeof(Scalar))/sizeof(ListEl)); 78 m_buffer = new Scalar[allocSize]; 79 } 80 else 81 { 82 m_allocatedElements = convert_index((size*sizeof(Scalar))/sizeof(ListEl)); 83 m_buffer = new Scalar[size]; 84 } 85 m_size = convert_index(size); 86 m_start = 0; 87 m_end = m_size; 88 } 89 90 void reallocateSparse() 91 { 92 Index copyElements = m_allocatedElements; 93 m_allocatedElements = (std::min)(StorageIndex(m_allocatedElements*1.5),m_size); 94 Index allocSize = m_allocatedElements * sizeof(ListEl); 95 allocSize = (allocSize + sizeof(Scalar) - 1)/sizeof(Scalar); 96 Scalar* newBuffer = new Scalar[allocSize]; 97 memcpy(newBuffer, m_buffer, copyElements * sizeof(ListEl)); 98 delete[] m_buffer; 99 m_buffer = newBuffer; 100 } 101 102 protected: 103 // element type of the linked list 104 struct ListEl 105 { 106 StorageIndex next; 107 StorageIndex index; 108 Scalar value; 109 }; 110 111 // used to store data in both mode 112 Scalar* m_buffer; 113 Scalar m_zero; 114 StorageIndex m_size; 115 StorageIndex m_start; 116 StorageIndex m_end; 117 StorageIndex m_allocatedSize; 118 StorageIndex m_allocatedElements; 119 StorageIndex m_mode; 120 121 // linked list mode 122 StorageIndex m_llStart; 123 StorageIndex m_llCurrent; 124 StorageIndex m_llSize; 125}; 126 127/** \returns the number of non zeros in the current sub vector */ 128template<typename _Scalar,typename _StorageIndex> 129Index AmbiVector<_Scalar,_StorageIndex>::nonZeros() const 130{ 131 if (m_mode==IsSparse) 132 return m_llSize; 133 else 134 return m_end - m_start; 135} 136 137template<typename _Scalar,typename _StorageIndex> 138void AmbiVector<_Scalar,_StorageIndex>::init(double estimatedDensity) 139{ 140 if (estimatedDensity>0.1) 141 init(IsDense); 142 else 143 init(IsSparse); 144} 145 146template<typename _Scalar,typename _StorageIndex> 147void AmbiVector<_Scalar,_StorageIndex>::init(int mode) 148{ 149 m_mode = mode; 150 if (m_mode==IsSparse) 151 { 152 m_llSize = 0; 153 m_llStart = -1; 154 } 155} 156 157/** Must be called whenever we might perform a write access 158 * with an index smaller than the previous one. 159 * 160 * Don't worry, this function is extremely cheap. 161 */ 162template<typename _Scalar,typename _StorageIndex> 163void AmbiVector<_Scalar,_StorageIndex>::restart() 164{ 165 m_llCurrent = m_llStart; 166} 167 168/** Set all coefficients of current subvector to zero */ 169template<typename _Scalar,typename _StorageIndex> 170void AmbiVector<_Scalar,_StorageIndex>::setZero() 171{ 172 if (m_mode==IsDense) 173 { 174 for (Index i=m_start; i<m_end; ++i) 175 m_buffer[i] = Scalar(0); 176 } 177 else 178 { 179 eigen_assert(m_mode==IsSparse); 180 m_llSize = 0; 181 m_llStart = -1; 182 } 183} 184 185template<typename _Scalar,typename _StorageIndex> 186_Scalar& AmbiVector<_Scalar,_StorageIndex>::coeffRef(Index i) 187{ 188 if (m_mode==IsDense) 189 return m_buffer[i]; 190 else 191 { 192 ListEl* EIGEN_RESTRICT llElements = reinterpret_cast<ListEl*>(m_buffer); 193 // TODO factorize the following code to reduce code generation 194 eigen_assert(m_mode==IsSparse); 195 if (m_llSize==0) 196 { 197 // this is the first element 198 m_llStart = 0; 199 m_llCurrent = 0; 200 ++m_llSize; 201 llElements[0].value = Scalar(0); 202 llElements[0].index = convert_index(i); 203 llElements[0].next = -1; 204 return llElements[0].value; 205 } 206 else if (i<llElements[m_llStart].index) 207 { 208 // this is going to be the new first element of the list 209 ListEl& el = llElements[m_llSize]; 210 el.value = Scalar(0); 211 el.index = convert_index(i); 212 el.next = m_llStart; 213 m_llStart = m_llSize; 214 ++m_llSize; 215 m_llCurrent = m_llStart; 216 return el.value; 217 } 218 else 219 { 220 StorageIndex nextel = llElements[m_llCurrent].next; 221 eigen_assert(i>=llElements[m_llCurrent].index && "you must call restart() before inserting an element with lower or equal index"); 222 while (nextel >= 0 && llElements[nextel].index<=i) 223 { 224 m_llCurrent = nextel; 225 nextel = llElements[nextel].next; 226 } 227 228 if (llElements[m_llCurrent].index==i) 229 { 230 // the coefficient already exists and we found it ! 231 return llElements[m_llCurrent].value; 232 } 233 else 234 { 235 if (m_llSize>=m_allocatedElements) 236 { 237 reallocateSparse(); 238 llElements = reinterpret_cast<ListEl*>(m_buffer); 239 } 240 eigen_internal_assert(m_llSize<m_allocatedElements && "internal error: overflow in sparse mode"); 241 // let's insert a new coefficient 242 ListEl& el = llElements[m_llSize]; 243 el.value = Scalar(0); 244 el.index = convert_index(i); 245 el.next = llElements[m_llCurrent].next; 246 llElements[m_llCurrent].next = m_llSize; 247 ++m_llSize; 248 return el.value; 249 } 250 } 251 } 252} 253 254template<typename _Scalar,typename _StorageIndex> 255_Scalar& AmbiVector<_Scalar,_StorageIndex>::coeff(Index i) 256{ 257 if (m_mode==IsDense) 258 return m_buffer[i]; 259 else 260 { 261 ListEl* EIGEN_RESTRICT llElements = reinterpret_cast<ListEl*>(m_buffer); 262 eigen_assert(m_mode==IsSparse); 263 if ((m_llSize==0) || (i<llElements[m_llStart].index)) 264 { 265 return m_zero; 266 } 267 else 268 { 269 Index elid = m_llStart; 270 while (elid >= 0 && llElements[elid].index<i) 271 elid = llElements[elid].next; 272 273 if (llElements[elid].index==i) 274 return llElements[m_llCurrent].value; 275 else 276 return m_zero; 277 } 278 } 279} 280 281/** Iterator over the nonzero coefficients */ 282template<typename _Scalar,typename _StorageIndex> 283class AmbiVector<_Scalar,_StorageIndex>::Iterator 284{ 285 public: 286 typedef _Scalar Scalar; 287 typedef typename NumTraits<Scalar>::Real RealScalar; 288 289 /** Default constructor 290 * \param vec the vector on which we iterate 291 * \param epsilon the minimal value used to prune zero coefficients. 292 * In practice, all coefficients having a magnitude smaller than \a epsilon 293 * are skipped. 294 */ 295 explicit Iterator(const AmbiVector& vec, const RealScalar& epsilon = 0) 296 : m_vector(vec) 297 { 298 using std::abs; 299 m_epsilon = epsilon; 300 m_isDense = m_vector.m_mode==IsDense; 301 if (m_isDense) 302 { 303 m_currentEl = 0; // this is to avoid a compilation warning 304 m_cachedValue = 0; // this is to avoid a compilation warning 305 m_cachedIndex = m_vector.m_start-1; 306 ++(*this); 307 } 308 else 309 { 310 ListEl* EIGEN_RESTRICT llElements = reinterpret_cast<ListEl*>(m_vector.m_buffer); 311 m_currentEl = m_vector.m_llStart; 312 while (m_currentEl>=0 && abs(llElements[m_currentEl].value)<=m_epsilon) 313 m_currentEl = llElements[m_currentEl].next; 314 if (m_currentEl<0) 315 { 316 m_cachedValue = 0; // this is to avoid a compilation warning 317 m_cachedIndex = -1; 318 } 319 else 320 { 321 m_cachedIndex = llElements[m_currentEl].index; 322 m_cachedValue = llElements[m_currentEl].value; 323 } 324 } 325 } 326 327 StorageIndex index() const { return m_cachedIndex; } 328 Scalar value() const { return m_cachedValue; } 329 330 operator bool() const { return m_cachedIndex>=0; } 331 332 Iterator& operator++() 333 { 334 using std::abs; 335 if (m_isDense) 336 { 337 do { 338 ++m_cachedIndex; 339 } while (m_cachedIndex<m_vector.m_end && abs(m_vector.m_buffer[m_cachedIndex])<=m_epsilon); 340 if (m_cachedIndex<m_vector.m_end) 341 m_cachedValue = m_vector.m_buffer[m_cachedIndex]; 342 else 343 m_cachedIndex=-1; 344 } 345 else 346 { 347 ListEl* EIGEN_RESTRICT llElements = reinterpret_cast<ListEl*>(m_vector.m_buffer); 348 do { 349 m_currentEl = llElements[m_currentEl].next; 350 } while (m_currentEl>=0 && abs(llElements[m_currentEl].value)<=m_epsilon); 351 if (m_currentEl<0) 352 { 353 m_cachedIndex = -1; 354 } 355 else 356 { 357 m_cachedIndex = llElements[m_currentEl].index; 358 m_cachedValue = llElements[m_currentEl].value; 359 } 360 } 361 return *this; 362 } 363 364 protected: 365 const AmbiVector& m_vector; // the target vector 366 StorageIndex m_currentEl; // the current element in sparse/linked-list mode 367 RealScalar m_epsilon; // epsilon used to prune zero coefficients 368 StorageIndex m_cachedIndex; // current coordinate 369 Scalar m_cachedValue; // current value 370 bool m_isDense; // mode of the vector 371}; 372 373} // end namespace internal 374 375} // end namespace Eigen 376 377#endif // EIGEN_AMBIVECTOR_H 378