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