Merge of itasc branch. Project files, scons and cmake should be working. Makefile...
[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     void resize(int rows, int cols)
263     {
264 //       std::cerr << this << " resize " << rows << "x" << cols << "\n";
265       const int outerSize = IsRowMajor ? rows : cols;
266       m_innerSize = IsRowMajor ? cols : rows;
267       m_data.clear();
268       if (m_outerSize != outerSize)
269       {
270         delete[] m_outerIndex;
271         m_outerIndex = new int [outerSize+1];
272         m_outerSize = outerSize;
273         memset(m_outerIndex, 0, (m_outerSize+1)*sizeof(int));
274       }
275     }
276     void resizeNonZeros(int size)
277     {
278       m_data.resize(size);
279     }
280
281     inline SparseMatrix()
282       : m_outerSize(-1), m_innerSize(0), m_outerIndex(0)
283     {
284       resize(0, 0);
285     }
286
287     inline SparseMatrix(int rows, int cols)
288       : m_outerSize(0), m_innerSize(0), m_outerIndex(0)
289     {
290       resize(rows, cols);
291     }
292
293     template<typename OtherDerived>
294     inline SparseMatrix(const SparseMatrixBase<OtherDerived>& other)
295       : m_outerSize(0), m_innerSize(0), m_outerIndex(0)
296     {
297       *this = other.derived();
298     }
299
300     inline SparseMatrix(const SparseMatrix& other)
301       : Base(), m_outerSize(0), m_innerSize(0), m_outerIndex(0)
302     {
303       *this = other.derived();
304     }
305
306     inline void swap(SparseMatrix& other)
307     {
308       //EIGEN_DBG_SPARSE(std::cout << "SparseMatrix:: swap\n");
309       std::swap(m_outerIndex, other.m_outerIndex);
310       std::swap(m_innerSize, other.m_innerSize);
311       std::swap(m_outerSize, other.m_outerSize);
312       m_data.swap(other.m_data);
313     }
314
315     inline SparseMatrix& operator=(const SparseMatrix& other)
316     {
317 //       std::cout << "SparseMatrix& operator=(const SparseMatrix& other)\n";
318       if (other.isRValue())
319       {
320         swap(other.const_cast_derived());
321       }
322       else
323       {
324         resize(other.rows(), other.cols());
325         memcpy(m_outerIndex, other.m_outerIndex, (m_outerSize+1)*sizeof(int));
326         m_data = other.m_data;
327       }
328       return *this;
329     }
330
331     template<typename OtherDerived>
332     inline SparseMatrix& operator=(const SparseMatrixBase<OtherDerived>& other)
333     {
334       const bool needToTranspose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit);
335       if (needToTranspose)
336       {
337         // two passes algorithm:
338         //  1 - compute the number of coeffs per dest inner vector
339         //  2 - do the actual copy/eval
340         // Since each coeff of the rhs has to be evaluated twice, let's evauluate it if needed
341         //typedef typename ei_nested<OtherDerived,2>::type OtherCopy;
342         typedef typename ei_eval<OtherDerived>::type OtherCopy;
343         typedef typename ei_cleantype<OtherCopy>::type _OtherCopy;
344         OtherCopy otherCopy(other.derived());
345
346         resize(other.rows(), other.cols());
347         Eigen::Map<VectorXi>(m_outerIndex,outerSize()).setZero();
348         // pass 1
349         // FIXME the above copy could be merged with that pass
350         for (int j=0; j<otherCopy.outerSize(); ++j)
351           for (typename _OtherCopy::InnerIterator it(otherCopy, j); it; ++it)
352             ++m_outerIndex[it.index()];
353
354         // prefix sum
355         int count = 0;
356         VectorXi positions(outerSize());
357         for (int j=0; j<outerSize(); ++j)
358         {
359           int tmp = m_outerIndex[j];
360           m_outerIndex[j] = count;
361           positions[j] = count;
362           count += tmp;
363         }
364         m_outerIndex[outerSize()] = count;
365         // alloc
366         m_data.resize(count);
367         // pass 2
368         for (int j=0; j<otherCopy.outerSize(); ++j)
369           for (typename _OtherCopy::InnerIterator it(otherCopy, j); it; ++it)
370           {
371             int pos = positions[it.index()]++;
372             m_data.index(pos) = j;
373             m_data.value(pos) = it.value();
374           }
375
376         return *this;
377       }
378       else
379       {
380         // there is no special optimization
381         return SparseMatrixBase<SparseMatrix>::operator=(other.derived());
382       }
383     }
384
385     friend std::ostream & operator << (std::ostream & s, const SparseMatrix& m)
386     {
387       EIGEN_DBG_SPARSE(
388         s << "Nonzero entries:\n";
389         for (int i=0; i<m.nonZeros(); ++i)
390         {
391           s << "(" << m.m_data.value(i) << "," << m.m_data.index(i) << ") ";
392         }
393         s << std::endl;
394         s << std::endl;
395         s << "Column pointers:\n";
396         for (int i=0; i<m.outerSize(); ++i)
397         {
398           s << m.m_outerIndex[i] << " ";
399         }
400         s << " $" << std::endl;
401         s << std::endl;
402       );
403       s << static_cast<const SparseMatrixBase<SparseMatrix>&>(m);
404       return s;
405     }
406
407     /** Destructor */
408     inline ~SparseMatrix()
409     {
410       delete[] m_outerIndex;
411     }
412 };
413
414 template<typename Scalar, int _Flags>
415 class SparseMatrix<Scalar,_Flags>::InnerIterator
416 {
417   public:
418     InnerIterator(const SparseMatrix& mat, int outer)
419       : m_matrix(mat), m_outer(outer), m_id(mat.m_outerIndex[outer]), m_start(m_id), m_end(mat.m_outerIndex[outer+1])
420     {}
421
422     template<unsigned int Added, unsigned int Removed>
423     InnerIterator(const Flagged<SparseMatrix,Added,Removed>& mat, int outer)
424       : m_matrix(mat._expression()), m_outer(outer), m_id(m_matrix.m_outerIndex[outer]),
425         m_start(m_id), m_end(m_matrix.m_outerIndex[outer+1])
426     {}
427
428     inline InnerIterator& operator++() { m_id++; return *this; }
429
430     inline Scalar value() const { return m_matrix.m_data.value(m_id); }
431     inline Scalar& valueRef() { return const_cast<Scalar&>(m_matrix.m_data.value(m_id)); }
432
433     inline int index() const { return m_matrix.m_data.index(m_id); }
434     inline int row() const { return IsRowMajor ? m_outer : index(); }
435     inline int col() const { return IsRowMajor ? index() : m_outer; }
436
437     inline operator bool() const { return (m_id < m_end) && (m_id>=m_start); }
438
439   protected:
440     const SparseMatrix& m_matrix;
441     const int m_outer;
442     int m_id;
443     const int m_start;
444     const int m_end;
445 };
446
447 #endif // EIGEN_SPARSEMATRIX_H