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