Merge of itasc branch. Project files, scons and cmake should be working. Makefile...
[blender.git] / extern / Eigen2 / Eigen / src / Core / Part.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 Benoit Jacob <jacob.benoit.1@gmail.com>
5 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
6 //
7 // Eigen is free software; you can redistribute it and/or
8 // modify it under the terms of the GNU Lesser General Public
9 // License as published by the Free Software Foundation; either
10 // version 3 of the License, or (at your option) any later version.
11 //
12 // Alternatively, you can redistribute it and/or
13 // modify it under the terms of the GNU General Public License as
14 // published by the Free Software Foundation; either version 2 of
15 // the License, or (at your option) any later version.
16 //
17 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
18 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
20 // GNU General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public
23 // License and a copy of the GNU General Public License along with
24 // Eigen. If not, see <http://www.gnu.org/licenses/>.
25
26 #ifndef EIGEN_PART_H
27 #define EIGEN_PART_H
28
29 /** \nonstableyet
30   * \class Part
31   *
32   * \brief Expression of a triangular matrix extracted from a given matrix
33   *
34   * \param MatrixType the type of the object in which we are taking the triangular part
35   * \param Mode the kind of triangular matrix expression to construct. Can be UpperTriangular, StrictlyUpperTriangular,
36   *             UnitUpperTriangular, LowerTriangular, StrictlyLowerTriangular, UnitLowerTriangular. This is in fact a bit field; it must have either
37   *             UpperTriangularBit or LowerTriangularBit, and additionnaly it may have either ZeroDiagBit or
38   *             UnitDiagBit.
39   *
40   * This class represents an expression of the upper or lower triangular part of
41   * a square matrix, possibly with a further assumption on the diagonal. It is the return type
42   * of MatrixBase::part() and most of the time this is the only way it is used.
43   *
44   * \sa MatrixBase::part()
45   */
46 template<typename MatrixType, unsigned int Mode>
47 struct ei_traits<Part<MatrixType, Mode> > : ei_traits<MatrixType>
48 {
49   typedef typename ei_nested<MatrixType>::type MatrixTypeNested;
50   typedef typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested;
51   enum {
52     Flags = (_MatrixTypeNested::Flags & (HereditaryBits) & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit))) | Mode,
53     CoeffReadCost = _MatrixTypeNested::CoeffReadCost
54   };
55 };
56
57 template<typename MatrixType, unsigned int Mode> class Part
58   : public MatrixBase<Part<MatrixType, Mode> >
59 {
60   public:
61
62     EIGEN_GENERIC_PUBLIC_INTERFACE(Part)
63
64     inline Part(const MatrixType& matrix) : m_matrix(matrix)
65     { ei_assert(ei_are_flags_consistent<Mode>::ret); }
66
67     /** \sa MatrixBase::operator+=() */
68     template<typename Other> Part&  operator+=(const Other& other);
69     /** \sa MatrixBase::operator-=() */
70     template<typename Other> Part&  operator-=(const Other& other);
71     /** \sa MatrixBase::operator*=() */
72     Part&  operator*=(const typename ei_traits<MatrixType>::Scalar& other);
73     /** \sa MatrixBase::operator/=() */
74     Part&  operator/=(const typename ei_traits<MatrixType>::Scalar& other);
75
76     /** \sa operator=(), MatrixBase::lazyAssign() */
77     template<typename Other> void lazyAssign(const Other& other);
78     /** \sa MatrixBase::operator=() */
79     template<typename Other> Part& operator=(const Other& other);
80
81     inline int rows() const { return m_matrix.rows(); }
82     inline int cols() const { return m_matrix.cols(); }
83     inline int stride() const { return m_matrix.stride(); }
84
85     inline Scalar coeff(int row, int col) const
86     {
87       // SelfAdjointBit doesn't play any role here: just because a matrix is selfadjoint doesn't say anything about
88       // each individual coefficient, except for the not-very-useful-here fact that diagonal coefficients are real.
89       if( ((Flags & LowerTriangularBit) && (col>row)) || ((Flags & UpperTriangularBit) && (row>col)) )
90         return (Scalar)0;
91       if(Flags & UnitDiagBit)
92         return col==row ? (Scalar)1 : m_matrix.coeff(row, col);
93       else if(Flags & ZeroDiagBit)
94         return col==row ? (Scalar)0 : m_matrix.coeff(row, col);
95       else
96         return m_matrix.coeff(row, col);
97     }
98
99     inline Scalar& coeffRef(int row, int col)
100     {
101       EIGEN_STATIC_ASSERT(!(Flags & UnitDiagBit), WRITING_TO_TRIANGULAR_PART_WITH_UNIT_DIAGONAL_IS_NOT_SUPPORTED)
102       EIGEN_STATIC_ASSERT(!(Flags & SelfAdjointBit), COEFFICIENT_WRITE_ACCESS_TO_SELFADJOINT_NOT_SUPPORTED)
103       ei_assert(   (Mode==UpperTriangular && col>=row)
104                 || (Mode==LowerTriangular && col<=row)
105                 || (Mode==StrictlyUpperTriangular && col>row)
106                 || (Mode==StrictlyLowerTriangular && col<row));
107       return m_matrix.const_cast_derived().coeffRef(row, col);
108     }
109
110     /** \internal */
111     const MatrixType& _expression() const { return m_matrix; }
112
113     /** discard any writes to a row */
114     const Block<Part, 1, ColsAtCompileTime> row(int i) { return Base::row(i); }
115     const Block<Part, 1, ColsAtCompileTime> row(int i) const { return Base::row(i); }
116     /** discard any writes to a column */
117     const Block<Part, RowsAtCompileTime, 1> col(int i) { return Base::col(i); }
118     const Block<Part, RowsAtCompileTime, 1> col(int i) const { return Base::col(i); }
119
120     template<typename OtherDerived>
121     void swap(const MatrixBase<OtherDerived>& other)
122     {
123       Part<SwapWrapper<MatrixType>,Mode>(const_cast<MatrixType&>(m_matrix)).lazyAssign(other.derived());
124     }
125
126   protected:
127
128     const typename MatrixType::Nested m_matrix;
129 };
130
131 /** \nonstableyet
132   * \returns an expression of a triangular matrix extracted from the current matrix
133   *
134   * The parameter \a Mode can have the following values: \c UpperTriangular, \c StrictlyUpperTriangular, \c UnitUpperTriangular,
135   * \c LowerTriangular, \c StrictlyLowerTriangular, \c UnitLowerTriangular.
136   *
137   * \addexample PartExample \label How to extract a triangular part of an arbitrary matrix
138   *
139   * Example: \include MatrixBase_extract.cpp
140   * Output: \verbinclude MatrixBase_extract.out
141   *
142   * \sa class Part, part(), marked()
143   */
144 template<typename Derived>
145 template<unsigned int Mode>
146 const Part<Derived, Mode> MatrixBase<Derived>::part() const
147 {
148   return derived();
149 }
150
151 template<typename MatrixType, unsigned int Mode>
152 template<typename Other>
153 inline Part<MatrixType, Mode>& Part<MatrixType, Mode>::operator=(const Other& other)
154 {
155   if(Other::Flags & EvalBeforeAssigningBit)
156   {
157     typename MatrixBase<Other>::PlainMatrixType other_evaluated(other.rows(), other.cols());
158     other_evaluated.template part<Mode>().lazyAssign(other);
159     lazyAssign(other_evaluated);
160   }
161   else
162     lazyAssign(other.derived());
163   return *this;
164 }
165
166 template<typename Derived1, typename Derived2, unsigned int Mode, int UnrollCount>
167 struct ei_part_assignment_impl
168 {
169   enum {
170     col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
171     row = (UnrollCount-1) % Derived1::RowsAtCompileTime
172   };
173
174   inline static void run(Derived1 &dst, const Derived2 &src)
175   {
176     ei_part_assignment_impl<Derived1, Derived2, Mode, UnrollCount-1>::run(dst, src);
177
178     if(Mode == SelfAdjoint)
179     {
180       if(row == col)
181         dst.coeffRef(row, col) = ei_real(src.coeff(row, col));
182       else if(row < col)
183         dst.coeffRef(col, row) = ei_conj(dst.coeffRef(row, col) = src.coeff(row, col));
184     }
185     else
186     {
187       ei_assert(Mode == UpperTriangular || Mode == LowerTriangular || Mode == StrictlyUpperTriangular || Mode == StrictlyLowerTriangular);
188       if((Mode == UpperTriangular && row <= col)
189       || (Mode == LowerTriangular && row >= col)
190       || (Mode == StrictlyUpperTriangular && row < col)
191       || (Mode == StrictlyLowerTriangular && row > col))
192         dst.copyCoeff(row, col, src);
193     }
194   }
195 };
196
197 template<typename Derived1, typename Derived2, unsigned int Mode>
198 struct ei_part_assignment_impl<Derived1, Derived2, Mode, 1>
199 {
200   inline static void run(Derived1 &dst, const Derived2 &src)
201   {
202     if(!(Mode & ZeroDiagBit))
203       dst.copyCoeff(0, 0, src);
204   }
205 };
206
207 // prevent buggy user code from causing an infinite recursion
208 template<typename Derived1, typename Derived2, unsigned int Mode>
209 struct ei_part_assignment_impl<Derived1, Derived2, Mode, 0>
210 {
211   inline static void run(Derived1 &, const Derived2 &) {}
212 };
213
214 template<typename Derived1, typename Derived2>
215 struct ei_part_assignment_impl<Derived1, Derived2, UpperTriangular, Dynamic>
216 {
217   inline static void run(Derived1 &dst, const Derived2 &src)
218   {
219     for(int j = 0; j < dst.cols(); ++j)
220       for(int i = 0; i <= j; ++i)
221         dst.copyCoeff(i, j, src);
222   }
223 };
224
225 template<typename Derived1, typename Derived2>
226 struct ei_part_assignment_impl<Derived1, Derived2, LowerTriangular, Dynamic>
227 {
228   inline static void run(Derived1 &dst, const Derived2 &src)
229   {
230     for(int j = 0; j < dst.cols(); ++j)
231       for(int i = j; i < dst.rows(); ++i)
232         dst.copyCoeff(i, j, src);
233   }
234 };
235
236 template<typename Derived1, typename Derived2>
237 struct ei_part_assignment_impl<Derived1, Derived2, StrictlyUpperTriangular, Dynamic>
238 {
239   inline static void run(Derived1 &dst, const Derived2 &src)
240   {
241     for(int j = 0; j < dst.cols(); ++j)
242       for(int i = 0; i < j; ++i)
243         dst.copyCoeff(i, j, src);
244   }
245 };
246 template<typename Derived1, typename Derived2>
247 struct ei_part_assignment_impl<Derived1, Derived2, StrictlyLowerTriangular, Dynamic>
248 {
249   inline static void run(Derived1 &dst, const Derived2 &src)
250   {
251     for(int j = 0; j < dst.cols(); ++j)
252       for(int i = j+1; i < dst.rows(); ++i)
253         dst.copyCoeff(i, j, src);
254   }
255 };
256 template<typename Derived1, typename Derived2>
257 struct ei_part_assignment_impl<Derived1, Derived2, SelfAdjoint, Dynamic>
258 {
259   inline static void run(Derived1 &dst, const Derived2 &src)
260   {
261     for(int j = 0; j < dst.cols(); ++j)
262     {
263       for(int i = 0; i < j; ++i)
264         dst.coeffRef(j, i) = ei_conj(dst.coeffRef(i, j) = src.coeff(i, j));
265       dst.coeffRef(j, j) = ei_real(src.coeff(j, j));
266     }
267   }
268 };
269
270 template<typename MatrixType, unsigned int Mode>
271 template<typename Other>
272 void Part<MatrixType, Mode>::lazyAssign(const Other& other)
273 {
274   const bool unroll = MatrixType::SizeAtCompileTime * Other::CoeffReadCost / 2 <= EIGEN_UNROLLING_LIMIT;
275   ei_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols());
276
277   ei_part_assignment_impl
278     <MatrixType, Other, Mode,
279     unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic
280     >::run(m_matrix.const_cast_derived(), other.derived());
281 }
282
283 /** \nonstableyet
284   * \returns a lvalue pseudo-expression allowing to perform special operations on \c *this.
285   *
286   * The \a Mode parameter can have the following values: \c UpperTriangular, \c StrictlyUpperTriangular, \c LowerTriangular,
287   * \c StrictlyLowerTriangular, \c SelfAdjoint.
288   *
289   * \addexample PartExample \label How to write to a triangular part of a matrix
290   *
291   * Example: \include MatrixBase_part.cpp
292   * Output: \verbinclude MatrixBase_part.out
293   *
294   * \sa class Part, MatrixBase::extract(), MatrixBase::marked()
295   */
296 template<typename Derived>
297 template<unsigned int Mode>
298 inline Part<Derived, Mode> MatrixBase<Derived>::part()
299 {
300   return Part<Derived, Mode>(derived());
301 }
302
303 /** \returns true if *this is approximately equal to an upper triangular matrix,
304   *          within the precision given by \a prec.
305   *
306   * \sa isLowerTriangular(), extract(), part(), marked()
307   */
308 template<typename Derived>
309 bool MatrixBase<Derived>::isUpperTriangular(RealScalar prec) const
310 {
311   if(cols() != rows()) return false;
312   RealScalar maxAbsOnUpperTriangularPart = static_cast<RealScalar>(-1);
313   for(int j = 0; j < cols(); ++j)
314     for(int i = 0; i <= j; ++i)
315     {
316       RealScalar absValue = ei_abs(coeff(i,j));
317       if(absValue > maxAbsOnUpperTriangularPart) maxAbsOnUpperTriangularPart = absValue;
318     }
319   for(int j = 0; j < cols()-1; ++j)
320     for(int i = j+1; i < rows(); ++i)
321       if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnUpperTriangularPart, prec)) return false;
322   return true;
323 }
324
325 /** \returns true if *this is approximately equal to a lower triangular matrix,
326   *          within the precision given by \a prec.
327   *
328   * \sa isUpperTriangular(), extract(), part(), marked()
329   */
330 template<typename Derived>
331 bool MatrixBase<Derived>::isLowerTriangular(RealScalar prec) const
332 {
333   if(cols() != rows()) return false;
334   RealScalar maxAbsOnLowerTriangularPart = static_cast<RealScalar>(-1);
335   for(int j = 0; j < cols(); ++j)
336     for(int i = j; i < rows(); ++i)
337     {
338       RealScalar absValue = ei_abs(coeff(i,j));
339       if(absValue > maxAbsOnLowerTriangularPart) maxAbsOnLowerTriangularPart = absValue;
340     }
341   for(int j = 1; j < cols(); ++j)
342     for(int i = 0; i < j; ++i)
343       if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnLowerTriangularPart, prec)) return false;
344   return true;
345 }
346
347 template<typename MatrixType, unsigned int Mode>
348 template<typename Other>
349 inline Part<MatrixType, Mode>& Part<MatrixType, Mode>::operator+=(const Other& other)
350 {
351   return *this = m_matrix + other;
352 }
353
354 template<typename MatrixType, unsigned int Mode>
355 template<typename Other>
356 inline Part<MatrixType, Mode>& Part<MatrixType, Mode>::operator-=(const Other& other)
357 {
358   return *this = m_matrix - other;
359 }
360
361 template<typename MatrixType, unsigned int Mode>
362 inline Part<MatrixType, Mode>& Part<MatrixType, Mode>::operator*=
363 (const typename ei_traits<MatrixType>::Scalar& other)
364 {
365   return *this = m_matrix * other;
366 }
367
368 template<typename MatrixType, unsigned int Mode>
369 inline Part<MatrixType, Mode>& Part<MatrixType, Mode>::operator/=
370 (const typename ei_traits<MatrixType>::Scalar& other)
371 {
372   return *this = m_matrix / other;
373 }
374
375 #endif // EIGEN_PART_H