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
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
AmbiVector(Index size)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 */
setBounds(Index start,Index end)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
~AmbiVector()52 ~AmbiVector() { delete[] m_buffer; }
53
resize(Index size)54 void resize(Index size)
55 {
56 if (m_allocatedSize < size)
57 reallocate(size);
58 m_size = size;
59 }
60
size()61 Index size() const { return m_size; }
62
63 protected:
64
reallocate(Index size)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) - 1)/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
reallocateSparse()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) - 1)/sizeof(Scalar);
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>
nonZeros()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>
init(double estimatedDensity)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>
init(int mode)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>
restart()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>
setZero()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>
coeffRef(_Index i)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>
coeff(_Index i)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)
m_vector(vec)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
index()323 Index index() const { return m_cachedIndex; }
value()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