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 std::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