Merge of itasc branch. Project files, scons and cmake should be working. Makefile...
[blender.git] / extern / Eigen2 / Eigen / src / Sparse / DynamicSparseMatrix.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_DYNAMIC_SPARSEMATRIX_H
26 #define EIGEN_DYNAMIC_SPARSEMATRIX_H
27
28 /** \class DynamicSparseMatrix
29   *
30   * \brief A sparse matrix class designed for matrix assembly purpose
31   *
32   * \param _Scalar the scalar type, i.e. the type of the coefficients
33   *
34   * Unlike SparseMatrix, this class provides a much higher degree of flexibility. In particular, it allows
35   * random read/write accesses in log(rho*outer_size) where \c rho is the probability that a coefficient is
36   * nonzero and outer_size is the number of columns if the matrix is column-major and the number of rows
37   * otherwise.
38   * 
39   * Internally, the data are stored as a std::vector of compressed vector. The performances of random writes might
40   * decrease as the number of nonzeros per inner-vector increase. In practice, we observed very good performance
41   * till about 100 nonzeros/vector, and the performance remains relatively good till 500 nonzeros/vectors.
42   *
43   * \see SparseMatrix
44   */
45 template<typename _Scalar, int _Flags>
46 struct ei_traits<DynamicSparseMatrix<_Scalar, _Flags> >
47 {
48   typedef _Scalar Scalar;
49   enum {
50     RowsAtCompileTime = Dynamic,
51     ColsAtCompileTime = Dynamic,
52     MaxRowsAtCompileTime = Dynamic,
53     MaxColsAtCompileTime = Dynamic,
54     Flags = SparseBit | _Flags,
55     CoeffReadCost = NumTraits<Scalar>::ReadCost,
56     SupportedAccessPatterns = OuterRandomAccessPattern
57   };
58 };
59
60 template<typename _Scalar, int _Flags>
61 class DynamicSparseMatrix
62   : public SparseMatrixBase<DynamicSparseMatrix<_Scalar, _Flags> >
63 {
64   public:
65     EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(DynamicSparseMatrix)
66     // FIXME: why are these operator already alvailable ???
67     // EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(DynamicSparseMatrix, +=)
68     // EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(DynamicSparseMatrix, -=)
69     typedef MappedSparseMatrix<Scalar,Flags> Map;
70
71   protected:
72
73     enum { IsRowMajor = Base::IsRowMajor };
74     typedef DynamicSparseMatrix<Scalar,(Flags&~RowMajorBit)|(IsRowMajor?RowMajorBit:0)> TransposedSparseMatrix;
75
76     int m_innerSize;
77     std::vector<CompressedStorage<Scalar> > m_data;
78
79   public:
80   
81     inline int rows() const { return IsRowMajor ? outerSize() : m_innerSize; }
82     inline int cols() const { return IsRowMajor ? m_innerSize : outerSize(); }
83     inline int innerSize() const { return m_innerSize; }
84     inline int outerSize() const { return m_data.size(); }
85     inline int innerNonZeros(int j) const { return m_data[j].size(); }
86     
87     std::vector<CompressedStorage<Scalar> >& _data() { return m_data; }
88     const std::vector<CompressedStorage<Scalar> >& _data() const { return m_data; }
89
90     /** \returns the coefficient value at given position \a row, \a col
91       * This operation involes a log(rho*outer_size) binary search.
92       */
93     inline Scalar coeff(int row, int col) const
94     {
95       const int outer = IsRowMajor ? row : col;
96       const int inner = IsRowMajor ? col : row;
97       return m_data[outer].at(inner);
98     }
99
100     /** \returns a reference to the coefficient value at given position \a row, \a col
101       * This operation involes a log(rho*outer_size) binary search. If the coefficient does not
102       * exist yet, then a sorted insertion into a sequential buffer is performed.
103       */
104     inline Scalar& coeffRef(int row, int col)
105     {
106       const int outer = IsRowMajor ? row : col;
107       const int inner = IsRowMajor ? col : row;
108       return m_data[outer].atWithInsertion(inner);
109     }
110
111     class InnerIterator;
112
113     inline void setZero()
114     {
115       for (int j=0; j<outerSize(); ++j)
116         m_data[j].clear();
117     }
118
119     /** \returns the number of non zero coefficients */
120     inline int nonZeros() const
121     {
122       int res = 0;
123       for (int j=0; j<outerSize(); ++j)
124         res += m_data[j].size();
125       return res;
126     }
127
128     /** Set the matrix to zero and reserve the memory for \a reserveSize nonzero coefficients. */
129     inline void startFill(int reserveSize = 1000)
130     {
131       if (outerSize()>0)
132       {
133         int reserveSizePerVector = std::max(reserveSize/outerSize(),4);
134         for (int j=0; j<outerSize(); ++j)
135         {
136           m_data[j].clear();
137           m_data[j].reserve(reserveSizePerVector);
138         }
139       }
140     }
141
142     /** inserts a nonzero coefficient at given coordinates \a row, \a col and returns its reference assuming that:
143       *  1 - the coefficient does not exist yet
144       *  2 - this the coefficient with greater inner coordinate for the given outer coordinate.
145       * In other words, assuming \c *this is column-major, then there must not exists any nonzero coefficient of coordinates
146       * \c i \c x \a col such that \c i >= \a row. Otherwise the matrix is invalid.
147       * 
148       * \see fillrand(), coeffRef()
149       */
150     inline Scalar& fill(int row, int col)
151     {
152       const int outer = IsRowMajor ? row : col;
153       const int inner = IsRowMajor ? col : row;
154       ei_assert(outer<int(m_data.size()) && inner<m_innerSize);
155       ei_assert((m_data[outer].size()==0) || (m_data[outer].index(m_data[outer].size()-1)<inner));
156       m_data[outer].append(0, inner);
157       return m_data[outer].value(m_data[outer].size()-1);
158     }
159
160     /** Like fill() but with random inner coordinates.
161       * Compared to the generic coeffRef(), the unique limitation is that we assume
162       * the coefficient does not exist yet.
163       */
164     inline Scalar& fillrand(int row, int col)
165     {
166       const int outer = IsRowMajor ? row : col;
167       const int inner = IsRowMajor ? col : row;
168       
169       int startId = 0;
170       int id = m_data[outer].size() - 1;
171       m_data[outer].resize(id+2,1);
172
173       while ( (id >= startId) && (m_data[outer].index(id) > inner) )
174       {
175         m_data[outer].index(id+1) = m_data[outer].index(id);
176         m_data[outer].value(id+1) = m_data[outer].value(id);
177         --id;
178       }
179       m_data[outer].index(id+1) = inner;
180       m_data[outer].value(id+1) = 0;
181       return m_data[outer].value(id+1);
182     }
183
184     /** Does nothing. Provided for compatibility with SparseMatrix. */
185     inline void endFill() {}
186     
187     void prune(Scalar reference, RealScalar epsilon = precision<RealScalar>())
188     {
189       for (int j=0; j<outerSize(); ++j)
190         m_data[j].prune(reference,epsilon);
191     }
192
193     /** Resize the matrix without preserving the data (the matrix is set to zero)
194       */
195     void resize(int rows, int cols)
196     {
197       const int outerSize = IsRowMajor ? rows : cols;
198       m_innerSize = IsRowMajor ? cols : rows;
199       setZero();
200       if (int(m_data.size()) != outerSize)
201       {
202         m_data.resize(outerSize);
203       }
204     }
205     
206     void resizeAndKeepData(int rows, int cols)
207     {
208       const int outerSize = IsRowMajor ? rows : cols;
209       const int innerSize = IsRowMajor ? cols : rows;
210       if (m_innerSize>innerSize)
211       {
212         // remove all coefficients with innerCoord>=innerSize
213         // TODO
214         std::cerr << "not implemented yet\n";
215         exit(2);
216       }
217       if (m_data.size() != outerSize)
218       {
219         m_data.resize(outerSize);
220       }
221     }
222
223     inline DynamicSparseMatrix()
224       : m_innerSize(0), m_data(0)
225     {
226       ei_assert(innerSize()==0 && outerSize()==0);
227     }
228
229     inline DynamicSparseMatrix(int rows, int cols)
230       : m_innerSize(0)
231     {
232       resize(rows, cols);
233     }
234
235     template<typename OtherDerived>
236     inline DynamicSparseMatrix(const SparseMatrixBase<OtherDerived>& other)
237       : m_innerSize(0)
238     {
239       *this = other.derived();
240     }
241
242     inline DynamicSparseMatrix(const DynamicSparseMatrix& other)
243       : Base(), m_innerSize(0)
244     {
245       *this = other.derived();
246     }
247
248     inline void swap(DynamicSparseMatrix& other)
249     {
250       //EIGEN_DBG_SPARSE(std::cout << "SparseMatrix:: swap\n");
251       std::swap(m_innerSize, other.m_innerSize);
252       //std::swap(m_outerSize, other.m_outerSize);
253       m_data.swap(other.m_data);
254     }
255
256     inline DynamicSparseMatrix& operator=(const DynamicSparseMatrix& other)
257     {
258       if (other.isRValue())
259       {
260         swap(other.const_cast_derived());
261       }
262       else
263       {
264         resize(other.rows(), other.cols());
265         m_data = other.m_data;
266       }
267       return *this;
268     }
269
270     template<typename OtherDerived>
271     inline DynamicSparseMatrix& operator=(const SparseMatrixBase<OtherDerived>& other)
272     {
273       return SparseMatrixBase<DynamicSparseMatrix>::operator=(other.derived());
274     }
275
276     /** Destructor */
277     inline ~DynamicSparseMatrix() {}
278 };
279
280 template<typename Scalar, int _Flags>
281 class DynamicSparseMatrix<Scalar,_Flags>::InnerIterator : public SparseVector<Scalar,_Flags>::InnerIterator
282 {
283     typedef typename SparseVector<Scalar,_Flags>::InnerIterator Base;
284   public:
285     InnerIterator(const DynamicSparseMatrix& mat, int outer)
286       : Base(mat.m_data[outer]), m_outer(outer)
287     {}
288     
289     inline int row() const { return IsRowMajor ? m_outer : Base::index(); }
290     inline int col() const { return IsRowMajor ? Base::index() : m_outer; }
291
292
293   protected:
294     const int m_outer;
295 };
296
297 #endif // EIGEN_DYNAMIC_SPARSEMATRIX_H