Merge of itasc branch. Project files, scons and cmake should be working. Makefile...
[blender.git] / extern / Eigen2 / Eigen / src / Sparse / SparseVector.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-2009 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_SPARSEVECTOR_H
26 #define EIGEN_SPARSEVECTOR_H
27
28 /** \class SparseVector
29   *
30   * \brief a sparse vector class
31   *
32   * \param _Scalar the scalar type, i.e. the type of the coefficients
33   *
34   * See http://www.netlib.org/linalg/html_templates/node91.html for details on the storage scheme.
35   *
36   */
37 template<typename _Scalar, int _Flags>
38 struct ei_traits<SparseVector<_Scalar, _Flags> >
39 {
40   typedef _Scalar Scalar;
41   enum {
42     IsColVector = _Flags & RowMajorBit ? 0 : 1,
43
44     RowsAtCompileTime = IsColVector ? Dynamic : 1,
45     ColsAtCompileTime = IsColVector ? 1 : Dynamic,
46     MaxRowsAtCompileTime = RowsAtCompileTime,
47     MaxColsAtCompileTime = ColsAtCompileTime,
48     Flags = SparseBit | _Flags,
49     CoeffReadCost = NumTraits<Scalar>::ReadCost,
50     SupportedAccessPatterns = InnerRandomAccessPattern
51   };
52 };
53
54 template<typename _Scalar, int _Flags>
55 class SparseVector
56   : public SparseMatrixBase<SparseVector<_Scalar, _Flags> >
57 {
58   public:
59     EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseVector)
60     EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(SparseVector, +=)
61     EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(SparseVector, -=)
62 //     EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(SparseVector, =)
63
64   protected:
65   public:
66
67     typedef SparseMatrixBase<SparseVector> SparseBase;
68     enum { IsColVector = ei_traits<SparseVector>::IsColVector };
69
70     CompressedStorage<Scalar> m_data;
71     int m_size;
72     
73     CompressedStorage<Scalar>& _data() { return m_data; }
74     CompressedStorage<Scalar>& _data() const { return m_data; }
75
76   public:
77
78     EIGEN_STRONG_INLINE int rows() const { return IsColVector ? m_size : 1; }
79     EIGEN_STRONG_INLINE int cols() const { return IsColVector ? 1 : m_size; }
80     EIGEN_STRONG_INLINE int innerSize() const { return m_size; }
81     EIGEN_STRONG_INLINE int outerSize() const { return 1; }
82     EIGEN_STRONG_INLINE int innerNonZeros(int j) const { ei_assert(j==0); return m_size; }
83
84     EIGEN_STRONG_INLINE const Scalar* _valuePtr() const { return &m_data.value(0); }
85     EIGEN_STRONG_INLINE Scalar* _valuePtr() { return &m_data.value(0); }
86
87     EIGEN_STRONG_INLINE const int* _innerIndexPtr() const { return &m_data.index(0); }
88     EIGEN_STRONG_INLINE int* _innerIndexPtr() { return &m_data.index(0); }
89
90     inline Scalar coeff(int row, int col) const
91     {
92       ei_assert((IsColVector ? col : row)==0);
93       return coeff(IsColVector ? row : col);
94     }
95     inline Scalar coeff(int i) const { return m_data.at(i); }
96
97     inline Scalar& coeffRef(int row, int col)
98     {
99       ei_assert((IsColVector ? col : row)==0);
100       return coeff(IsColVector ? row : col);
101     }
102
103     /** \returns a reference to the coefficient value at given index \a i
104       * This operation involes a log(rho*size) binary search. If the coefficient does not
105       * exist yet, then a sorted insertion into a sequential buffer is performed.
106       * 
107       * This insertion might be very costly if the number of nonzeros above \a i is large.
108       */
109     inline Scalar& coeffRef(int i)
110     {
111       return m_data.atWithInsertion(i);
112     }
113
114   public:
115
116     class InnerIterator;
117
118     inline void setZero() { m_data.clear(); }
119
120     /** \returns the number of non zero coefficients */
121     inline int nonZeros() const  { return m_data.size(); }
122
123     /**
124       */
125     inline void reserve(int reserveSize) { m_data.reserve(reserveSize); }
126     
127     inline void startFill(int reserve)
128     {
129       setZero();
130       m_data.reserve(reserve);
131     }
132
133     /**
134       */
135     inline Scalar& fill(int r, int c)
136     {
137       ei_assert(r==0 || c==0);
138       return fill(IsColVector ? r : c);
139     }
140     
141     inline Scalar& fill(int i)
142     {
143       m_data.append(0, i);
144       return m_data.value(m_data.size()-1);
145     }
146     
147     inline Scalar& fillrand(int r, int c)
148     {
149       ei_assert(r==0 || c==0);
150       return fillrand(IsColVector ? r : c);
151     }
152
153     /** Like fill() but with random coordinates.
154       */
155     inline Scalar& fillrand(int i)
156     {
157       int startId = 0;
158       int id = m_data.size() - 1;
159       m_data.resize(id+2,1);
160
161       while ( (id >= startId) && (m_data.index(id) > i) )
162       {
163         m_data.index(id+1) = m_data.index(id);
164         m_data.value(id+1) = m_data.value(id);
165         --id;
166       }
167       m_data.index(id+1) = i;
168       m_data.value(id+1) = 0;
169       return m_data.value(id+1);
170     }
171     
172     inline void endFill() {}
173     
174     void prune(Scalar reference, RealScalar epsilon = precision<RealScalar>())
175     {
176       m_data.prune(reference,epsilon);
177     }
178     
179     void resize(int rows, int cols)
180     {
181       ei_assert(rows==1 || cols==1);
182       resize(IsColVector ? rows : cols);
183     }
184
185     void resize(int newSize)
186     {
187       m_size = newSize;
188       m_data.clear();
189     }
190
191     void resizeNonZeros(int size) { m_data.resize(size); }
192
193     inline SparseVector() : m_size(0) { resize(0); }
194
195     inline SparseVector(int size) : m_size(0) { resize(size); }
196     
197     inline SparseVector(int rows, int cols) : m_size(0) { resize(rows,cols); }
198
199     template<typename OtherDerived>
200     inline SparseVector(const MatrixBase<OtherDerived>& other)
201       : m_size(0)
202     {
203       *this = other.derived();
204     }
205     
206     template<typename OtherDerived>
207     inline SparseVector(const SparseMatrixBase<OtherDerived>& other)
208       : m_size(0)
209     {
210       *this = other.derived();
211     }
212
213     inline SparseVector(const SparseVector& other)
214       : m_size(0)
215     {
216       *this = other.derived();
217     }
218
219     inline void swap(SparseVector& other)
220     {
221       std::swap(m_size, other.m_size);
222       m_data.swap(other.m_data);
223     }
224
225     inline SparseVector& operator=(const SparseVector& other)
226     {
227       if (other.isRValue())
228       {
229         swap(other.const_cast_derived());
230       }
231       else
232       {
233         resize(other.size());
234         m_data = other.m_data;
235       }
236       return *this;
237     }
238
239     template<typename OtherDerived>
240     inline SparseVector& operator=(const SparseMatrixBase<OtherDerived>& other)
241     {
242       return Base::operator=(other);
243     }
244     
245 //       const bool needToTranspose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit);
246 //       if (needToTranspose)
247 //       {
248 //         // two passes algorithm:
249 //         //  1 - compute the number of coeffs per dest inner vector
250 //         //  2 - do the actual copy/eval
251 //         // Since each coeff of the rhs has to be evaluated twice, let's evauluate it if needed
252 //         typedef typename ei_nested<OtherDerived,2>::type OtherCopy;
253 //         OtherCopy otherCopy(other.derived());
254 //         typedef typename ei_cleantype<OtherCopy>::type _OtherCopy;
255 //
256 //         resize(other.rows(), other.cols());
257 //         Eigen::Map<VectorXi>(m_outerIndex,outerSize()).setZero();
258 //         // pass 1
259 //         // FIXME the above copy could be merged with that pass
260 //         for (int j=0; j<otherCopy.outerSize(); ++j)
261 //           for (typename _OtherCopy::InnerIterator it(otherCopy, j); it; ++it)
262 //             ++m_outerIndex[it.index()];
263 //
264 //         // prefix sum
265 //         int count = 0;
266 //         VectorXi positions(outerSize());
267 //         for (int j=0; j<outerSize(); ++j)
268 //         {
269 //           int tmp = m_outerIndex[j];
270 //           m_outerIndex[j] = count;
271 //           positions[j] = count;
272 //           count += tmp;
273 //         }
274 //         m_outerIndex[outerSize()] = count;
275 //         // alloc
276 //         m_data.resize(count);
277 //         // pass 2
278 //         for (int j=0; j<otherCopy.outerSize(); ++j)
279 //           for (typename _OtherCopy::InnerIterator it(otherCopy, j); it; ++it)
280 //           {
281 //             int pos = positions[it.index()]++;
282 //             m_data.index(pos) = j;
283 //             m_data.value(pos) = it.value();
284 //           }
285 //
286 //         return *this;
287 //       }
288 //       else
289 //       {
290 //         // there is no special optimization
291 //         return SparseMatrixBase<SparseMatrix>::operator=(other.derived());
292 //       }
293 //     }
294
295     friend std::ostream & operator << (std::ostream & s, const SparseVector& m)
296     {
297       for (unsigned int i=0; i<m.nonZeros(); ++i)
298         s << "(" << m.m_data.value(i) << "," << m.m_data.index(i) << ") ";
299       s << std::endl;
300       return s;
301     }
302
303     // this specialized version does not seems to be faster
304 //     Scalar dot(const SparseVector& other) const
305 //     {
306 //       int i=0, j=0;
307 //       Scalar res = 0;
308 //       asm("#begindot");
309 //       while (i<nonZeros() && j<other.nonZeros())
310 //       {
311 //         if (m_data.index(i)==other.m_data.index(j))
312 //         {
313 //           res += m_data.value(i) * ei_conj(other.m_data.value(j));
314 //           ++i; ++j;
315 //         }
316 //         else if (m_data.index(i)<other.m_data.index(j))
317 //           ++i;
318 //         else
319 //           ++j;
320 //       }
321 //       asm("#enddot");
322 //       return res;
323 //     }
324
325     /** Destructor */
326     inline ~SparseVector() {}
327 };
328
329 template<typename Scalar, int _Flags>
330 class SparseVector<Scalar,_Flags>::InnerIterator
331 {
332   public:
333     InnerIterator(const SparseVector& vec, int outer=0)
334       : m_data(vec.m_data), m_id(0), m_end(m_data.size())
335     {
336       ei_assert(outer==0);
337     }
338     
339     InnerIterator(const CompressedStorage<Scalar>& data)
340       : m_data(data), m_id(0), m_end(m_data.size())
341     {}
342
343     template<unsigned int Added, unsigned int Removed>
344     InnerIterator(const Flagged<SparseVector,Added,Removed>& vec, int outer)
345       : m_data(vec._expression().m_data), m_id(0), m_end(m_data.size())
346     {}
347
348     inline InnerIterator& operator++() { m_id++; return *this; }
349
350     inline Scalar value() const { return m_data.value(m_id); }
351     inline Scalar& valueRef() { return const_cast<Scalar&>(m_data.value(m_id)); }
352
353     inline int index() const { return m_data.index(m_id); }
354     inline int row() const { return IsColVector ? index() : 0; }
355     inline int col() const { return IsColVector ? 0 : index(); }
356
357     inline operator bool() const { return (m_id < m_end); }
358
359   protected:
360     const CompressedStorage<Scalar>& m_data;
361     int m_id;
362     const int m_end;
363 };
364
365 #endif // EIGEN_SPARSEVECTOR_H