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