AmbiVector.h revision 7faaa9f3f0df9d23790277834d426c3d992ac3ba
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 _Index>
23class 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
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 */
124template<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
133template<typename _Scalar,typename _Index>
134void AmbiVector<_Scalar,_Index>::init(double estimatedDensity)
135{
136  if (estimatedDensity>0.1)
137    init(IsDense);
138  else
139    init(IsSparse);
140}
141
142template<typename _Scalar,typename _Index>
143void 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  */
158template<typename _Scalar,typename _Index>
159void AmbiVector<_Scalar,_Index>::restart()
160{
161  m_llCurrent = m_llStart;
162}
163
164/** Set all coefficients of current subvector to zero */
165template<typename _Scalar,typename _Index>
166void 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
181template<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
250template<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 */
278template<typename _Scalar,typename _Index>
279class 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