Cycles: svn merge -r41225:41232 ^/trunk/blender
[blender.git] / extern / Eigen2 / Eigen / src / Sparse / SparseMatrix.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_SPARSEMATRIX_H
26 #define EIGEN_SPARSEMATRIX_H
27
28 /** \class SparseMatrix
29   *
30   * \brief Sparse matrix
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<SparseMatrix<_Scalar, _Flags> >
39 {
40   typedef _Scalar Scalar;
41   enum {
42     RowsAtCompileTime = Dynamic,
43     ColsAtCompileTime = Dynamic,
44     MaxRowsAtCompileTime = Dynamic,
45     MaxColsAtCompileTime = Dynamic,
46     Flags = SparseBit | _Flags,
47     CoeffReadCost = NumTraits<Scalar>::ReadCost,
48     SupportedAccessPatterns = InnerRandomAccessPattern
49   };
50 };
51
52
53
54 template<typename _Scalar, int _Flags>
55 class SparseMatrix
56   : public SparseMatrixBase<SparseMatrix<_Scalar, _Flags> >
57 {
58   public:
59     EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseMatrix)
60     EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(SparseMatrix, +=)
61     EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(SparseMatrix, -=)
62     // FIXME: why are these operator already alvailable ???
63     // EIGEN_SPARSE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(SparseMatrix, *=)
64     // EIGEN_SPARSE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(SparseMatrix, /=)
65
66     typedef MappedSparseMatrix<Scalar,Flags> Map;
67
68   protected:
69
70     enum { IsRowMajor = Base::IsRowMajor };
71     typedef SparseMatrix<Scalar,(Flags&~RowMajorBit)|(IsRowMajor?RowMajorBit:0)> TransposedSparseMatrix;
72
73     int m_outerSize;
74     int m_innerSize;
75     int* m_outerIndex;
76     CompressedStorage<Scalar> m_data;
77
78   public:
79
80     inline int rows() const { return IsRowMajor ? m_outerSize : m_innerSize; }
81     inline int cols() const { return IsRowMajor ? m_innerSize : m_outerSize; }
82
83     inline int innerSize() const { return m_innerSize; }
84     inline int outerSize() const { return m_outerSize; }
85     inline int innerNonZeros(int j) const { return m_outerIndex[j+1]-m_outerIndex[j]; }
86
87     inline const Scalar* _valuePtr() const { return &m_data.value(0); }
88     inline Scalar* _valuePtr() { return &m_data.value(0); }
89
90     inline const int* _innerIndexPtr() const { return &m_data.index(0); }
91     inline int* _innerIndexPtr() { return &m_data.index(0); }
92
93     inline const int* _outerIndexPtr() const { return m_outerIndex; }
94     inline int* _outerIndexPtr() { return m_outerIndex; }
95
96     inline Scalar coeff(int row, int col) const
97     {
98       const int outer = IsRowMajor ? row : col;
99       const int inner = IsRowMajor ? col : row;
100       return m_data.atInRange(m_outerIndex[outer], m_outerIndex[outer+1], inner);
101     }
102
103     inline Scalar& coeffRef(int row, int col)
104     {
105       const int outer = IsRowMajor ? row : col;
106       const int inner = IsRowMajor ? col : row;
107
108       int start = m_outerIndex[outer];
109       int end = m_outerIndex[outer+1];
110       ei_assert(end>=start && "you probably called coeffRef on a non finalized matrix");
111       ei_assert(end>start && "coeffRef cannot be called on a zero coefficient");
112       const int id = m_data.searchLowerIndex(start,end-1,inner);
113       ei_assert((id<end) && (m_data.index(id)==inner) && "coeffRef cannot be called on a zero coefficient");
114       return m_data.value(id);
115     }
116
117   public:
118
119     class InnerIterator;
120
121     inline void setZero()
122     {
123       m_data.clear();
124       //if (m_outerSize)
125       memset(m_outerIndex, 0, (m_outerSize+1)*sizeof(int));
126 //       for (int i=0; i<m_outerSize; ++i)
127 //         m_outerIndex[i] = 0;
128 //       if (m_outerSize)
129 //         m_outerIndex[i] = 0;
130     }
131
132     /** \returns the number of non zero coefficients */
133     inline int nonZeros() const  { return m_data.size(); }
134
135     /** Initializes the filling process of \c *this.
136       * \param reserveSize approximate number of nonzeros
137       * Note that the matrix \c *this is zero-ed.
138       */
139     inline void startFill(int reserveSize = 1000)
140     {
141       setZero();
142       m_data.reserve(reserveSize);
143     }
144
145     /**
146       */
147     inline Scalar& fill(int row, int col)
148     {
149       const int outer = IsRowMajor ? row : col;
150       const int inner = IsRowMajor ? col : row;
151
152       if (m_outerIndex[outer+1]==0)
153       {
154         // we start a new inner vector
155         int i = outer;
156         while (i>=0 && m_outerIndex[i]==0)
157         {
158           m_outerIndex[i] = m_data.size();
159           --i;
160         }
161         m_outerIndex[outer+1] = m_outerIndex[outer];
162       }
163       else
164       {
165         ei_assert(m_data.index(m_data.size()-1)<inner && "wrong sorted insertion");
166       }
167       assert(size_t(m_outerIndex[outer+1]) == m_data.size());
168       int id = m_outerIndex[outer+1];
169       ++m_outerIndex[outer+1];
170
171       m_data.append(0, inner);
172       return m_data.value(id);
173     }
174
175     /** Like fill() but with random inner coordinates.
176       */
177     inline Scalar& fillrand(int row, int col)
178     {
179       const int outer = IsRowMajor ? row : col;
180       const int inner = IsRowMajor ? col : row;
181       if (m_outerIndex[outer+1]==0)
182       {
183         // we start a new inner vector
184         // nothing special to do here
185         int i = outer;
186         while (i>=0 && m_outerIndex[i]==0)
187         {
188           m_outerIndex[i] = m_data.size();
189           --i;
190         }
191         m_outerIndex[outer+1] = m_outerIndex[outer];
192       }
193       assert(size_t(m_outerIndex[outer+1]) == m_data.size() && "invalid outer index");
194       size_t startId = m_outerIndex[outer];
195       // FIXME let's make sure sizeof(long int) == sizeof(size_t)
196       size_t id = m_outerIndex[outer+1];
197       ++m_outerIndex[outer+1];
198
199       float reallocRatio = 1;
200       if (m_data.allocatedSize()<id+1)
201       {
202         // we need to reallocate the data, to reduce multiple reallocations
203         // we use a smart resize algorithm based on the current filling ratio
204         // we use float to avoid overflows
205         float nnzEstimate = float(m_outerIndex[outer])*float(m_outerSize)/float(outer);
206         reallocRatio = (nnzEstimate-float(m_data.size()))/float(m_data.size());
207         // let's bounds the realloc ratio to
208         //   1) reduce multiple minor realloc when the matrix is almost filled
209         //   2) avoid to allocate too much memory when the matrix is almost empty
210         reallocRatio = std::min(std::max(reallocRatio,1.5f),8.f);
211       }
212       m_data.resize(id+1,reallocRatio);
213
214       while ( (id > startId) && (m_data.index(id-1) > inner) )
215       {
216         m_data.index(id) = m_data.index(id-1);
217         m_data.value(id) = m_data.value(id-1);
218         --id;
219       }
220
221       m_data.index(id) = inner;
222       return (m_data.value(id) = 0);
223     }
224
225     inline void endFill()
226     {
227       int size = m_data.size();
228       int i = m_outerSize;
229       // find the last filled column
230       while (i>=0 && m_outerIndex[i]==0)
231         --i;
232       ++i;
233       while (i<=m_outerSize)
234       {
235         m_outerIndex[i] = size;
236         ++i;
237       }
238     }
239
240     void prune(Scalar reference, RealScalar epsilon = precision<RealScalar>())
241     {
242       int k = 0;
243       for (int j=0; j<m_outerSize; ++j)
244       {
245         int previousStart = m_outerIndex[j];
246         m_outerIndex[j] = k;
247         int end = m_outerIndex[j+1];
248         for (int i=previousStart; i<end; ++i)
249         {
250           if (!ei_isMuchSmallerThan(m_data.value(i), reference, epsilon))
251           {
252             m_data.value(k) = m_data.value(i);
253             m_data.index(k) = m_data.index(i);
254             ++k;
255           }
256         }
257       }
258       m_outerIndex[m_outerSize] = k;
259       m_data.resize(k,0);
260     }
261
262     /** Resizes the matrix to a \a rows x \a cols matrix and initializes it to zero
263       * \sa resizeNonZeros(int), reserve(), setZero()
264       */
265     void resize(int rows, int cols)
266     {
267       const int outerSize = IsRowMajor ? rows : cols;
268       m_innerSize = IsRowMajor ? cols : rows;
269       m_data.clear();
270       if (m_outerSize != outerSize || m_outerSize==0)
271       {
272         delete[] m_outerIndex;
273         m_outerIndex = new int [outerSize+1];
274         m_outerSize = outerSize;
275       }
276       memset(m_outerIndex, 0, (m_outerSize+1)*sizeof(int));
277     }
278     void resizeNonZeros(int size)
279     {
280       m_data.resize(size);
281     }
282
283     inline SparseMatrix()
284       : m_outerSize(-1), m_innerSize(0), m_outerIndex(0)
285     {
286       resize(0, 0);
287     }
288
289     inline SparseMatrix(int rows, int cols)
290       : m_outerSize(0), m_innerSize(0), m_outerIndex(0)
291     {
292       resize(rows, cols);
293     }
294
295     template<typename OtherDerived>
296     inline SparseMatrix(const SparseMatrixBase<OtherDerived>& other)
297       : m_outerSize(0), m_innerSize(0), m_outerIndex(0)
298     {
299       *this = other.derived();
300     }
301
302     inline SparseMatrix(const SparseMatrix& other)
303       : Base(), m_outerSize(0), m_innerSize(0), m_outerIndex(0)
304     {
305       *this = other.derived();
306     }
307
308     inline void swap(SparseMatrix& other)
309     {
310       //EIGEN_DBG_SPARSE(std::cout << "SparseMatrix:: swap\n");
311       std::swap(m_outerIndex, other.m_outerIndex);
312       std::swap(m_innerSize, other.m_innerSize);
313       std::swap(m_outerSize, other.m_outerSize);
314       m_data.swap(other.m_data);
315     }
316
317     inline SparseMatrix& operator=(const SparseMatrix& other)
318     {
319 //       std::cout << "SparseMatrix& operator=(const SparseMatrix& other)\n";
320       if (other.isRValue())
321       {
322         swap(other.const_cast_derived());
323       }
324       else
325       {
326         resize(other.rows(), other.cols());
327         memcpy(m_outerIndex, other.m_outerIndex, (m_outerSize+1)*sizeof(int));
328         m_data = other.m_data;
329       }
330       return *this;
331     }
332
333     template<typename OtherDerived>
334     inline SparseMatrix& operator=(const SparseMatrixBase<OtherDerived>& other)
335     {
336       const bool needToTranspose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit);
337       if (needToTranspose)
338       {
339         // two passes algorithm:
340         //  1 - compute the number of coeffs per dest inner vector
341         //  2 - do the actual copy/eval
342         // Since each coeff of the rhs has to be evaluated twice, let's evauluate it if needed
343         //typedef typename ei_nested<OtherDerived,2>::type OtherCopy;
344         typedef typename ei_eval<OtherDerived>::type OtherCopy;
345         typedef typename ei_cleantype<OtherCopy>::type _OtherCopy;
346         OtherCopy otherCopy(other.derived());
347
348         resize(other.rows(), other.cols());
349         Eigen::Map<VectorXi>(m_outerIndex,outerSize()).setZero();
350         // pass 1
351         // FIXME the above copy could be merged with that pass
352         for (int j=0; j<otherCopy.outerSize(); ++j)
353           for (typename _OtherCopy::InnerIterator it(otherCopy, j); it; ++it)
354             ++m_outerIndex[it.index()];
355
356         // prefix sum
357         int count = 0;
358         VectorXi positions(outerSize());
359         for (int j=0; j<outerSize(); ++j)
360         {
361           int tmp = m_outerIndex[j];
362           m_outerIndex[j] = count;
363           positions[j] = count;
364           count += tmp;
365         }
366         m_outerIndex[outerSize()] = count;
367         // alloc
368         m_data.resize(count);
369         // pass 2
370         for (int j=0; j<otherCopy.outerSize(); ++j)
371           for (typename _OtherCopy::InnerIterator it(otherCopy, j); it; ++it)
372           {
373             int pos = positions[it.index()]++;
374             m_data.index(pos) = j;
375             m_data.value(pos) = it.value();
376           }
377
378         return *this;
379       }
380       else
381       {
382         // there is no special optimization
383         return SparseMatrixBase<SparseMatrix>::operator=(other.derived());
384       }
385     }
386
387     friend std::ostream & operator << (std::ostream & s, const SparseMatrix& m)
388     {
389       EIGEN_DBG_SPARSE(
390         s << "Nonzero entries:\n";
391         for (int i=0; i<m.nonZeros(); ++i)
392         {
393           s << "(" << m.m_data.value(i) << "," << m.m_data.index(i) << ") ";
394         }
395         s << std::endl;
396         s << std::endl;
397         s << "Column pointers:\n";
398         for (int i=0; i<m.outerSize(); ++i)
399         {
400           s << m.m_outerIndex[i] << " ";
401         }
402         s << " $" << std::endl;
403         s << std::endl;
404       );
405       s << static_cast<const SparseMatrixBase<SparseMatrix>&>(m);
406       return s;
407     }
408
409     /** Destructor */
410     inline ~SparseMatrix()
411     {
412       delete[] m_outerIndex;
413     }
414 };
415
416 template<typename Scalar, int _Flags>
417 class SparseMatrix<Scalar,_Flags>::InnerIterator
418 {
419   public:
420     InnerIterator(const SparseMatrix& mat, int outer)
421       : m_matrix(mat), m_outer(outer), m_id(mat.m_outerIndex[outer]), m_start(m_id), m_end(mat.m_outerIndex[outer+1])
422     {}
423
424     template<unsigned int Added, unsigned int Removed>
425     InnerIterator(const Flagged<SparseMatrix,Added,Removed>& mat, int outer)
426       : m_matrix(mat._expression()), m_outer(outer), m_id(m_matrix.m_outerIndex[outer]),
427         m_start(m_id), m_end(m_matrix.m_outerIndex[outer+1])
428     {}
429
430     inline InnerIterator& operator++() { m_id++; return *this; }
431
432     inline Scalar value() const { return m_matrix.m_data.value(m_id); }
433     inline Scalar& valueRef() { return const_cast<Scalar&>(m_matrix.m_data.value(m_id)); }
434
435     inline int index() const { return m_matrix.m_data.index(m_id); }
436     inline int row() const { return IsRowMajor ? m_outer : index(); }
437     inline int col() const { return IsRowMajor ? index() : m_outer; }
438
439     inline operator bool() const { return (m_id < m_end) && (m_id>=m_start); }
440
441   protected:
442     const SparseMatrix& m_matrix;
443     const int m_outer;
444     int m_id;
445     const int m_start;
446     const int m_end;
447
448   private:
449     InnerIterator& operator=(const InnerIterator&);
450 };
451
452 #endif // EIGEN_SPARSEMATRIX_H