← Back to branch summary

~siretart/ubuntu/utopic/blender/libav10

« back to all changes in this revision

Viewing changes to extern/Eigen2/Eigen/src/Core/Part.h

  • Committer: Package Import Robot
  • Author(s): Matteo F. Vescovi
  • Date: 2012-07-23 08:54:18 UTC
  • mfrom: (14.2.16 sid)
  • mto: (14.2.19 sid)
  • mto: This revision was merged to the branch mainline in revision 42.
  • Revision ID: package-import@ubuntu.com-20120723085418-9foz30v6afaf5ffs
Tags: 2.63a-2
http://bugs.debian.org/658075
* debian/: Cycles support added (Closes: #658075)
  For now, this top feature has been enabled only
  on [any-amd64 any-i386] architectures because
  of OpenImageIO failing on all others
* debian/: scripts installation path changed
  from /usr/lib to /usr/share:
  + debian/patches/: patchset re-worked for path changing
  + debian/control: "Breaks" field added on yafaray-exporter

Show diffs side-by-side

added added

removed removed

Lines of Context:
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
 
    const typename MatrixType::Nested m_matrix;
128
 
 
129
 
  private:
130
 
    Part& operator=(const Part&);
131
 
};
132
 
 
133
 
/** \nonstableyet
134
 
  * \returns an expression of a triangular matrix extracted from the current matrix
135
 
  *
136
 
  * The parameter \a Mode can have the following values: \c UpperTriangular, \c StrictlyUpperTriangular, \c UnitUpperTriangular,
137
 
  * \c LowerTriangular, \c StrictlyLowerTriangular, \c UnitLowerTriangular.
138
 
  *
139
 
  * \addexample PartExample \label How to extract a triangular part of an arbitrary matrix
140
 
  *
141
 
  * Example: \include MatrixBase_extract.cpp
142
 
  * Output: \verbinclude MatrixBase_extract.out
143
 
  *
144
 
  * \sa class Part, part(), marked()
145
 
  */
146
 
template<typename Derived>
147
 
template<unsigned int Mode>
148
 
const Part<Derived, Mode> MatrixBase<Derived>::part() const
149
 
{
150
 
  return derived();
151
 
}
152
 
 
153
 
template<typename MatrixType, unsigned int Mode>
154
 
template<typename Other>
155
 
inline Part<MatrixType, Mode>& Part<MatrixType, Mode>::operator=(const Other& other)
156
 
{
157
 
  if(Other::Flags & EvalBeforeAssigningBit)
158
 
  {
159
 
    typename MatrixBase<Other>::PlainMatrixType other_evaluated(other.rows(), other.cols());
160
 
    other_evaluated.template part<Mode>().lazyAssign(other);
161
 
    lazyAssign(other_evaluated);
162
 
  }
163
 
  else
164
 
    lazyAssign(other.derived());
165
 
  return *this;
166
 
}
167
 
 
168
 
template<typename Derived1, typename Derived2, unsigned int Mode, int UnrollCount>
169
 
struct ei_part_assignment_impl
170
 
{
171
 
  enum {
172
 
    col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
173
 
    row = (UnrollCount-1) % Derived1::RowsAtCompileTime
174
 
  };
175
 
 
176
 
  inline static void run(Derived1 &dst, const Derived2 &src)
177
 
  {
178
 
    ei_part_assignment_impl<Derived1, Derived2, Mode, UnrollCount-1>::run(dst, src);
179
 
 
180
 
    if(Mode == SelfAdjoint)
181
 
    {
182
 
      if(row == col)
183
 
        dst.coeffRef(row, col) = ei_real(src.coeff(row, col));
184
 
      else if(row < col)
185
 
        dst.coeffRef(col, row) = ei_conj(dst.coeffRef(row, col) = src.coeff(row, col));
186
 
    }
187
 
    else
188
 
    {
189
 
      ei_assert(Mode == UpperTriangular || Mode == LowerTriangular || Mode == StrictlyUpperTriangular || Mode == StrictlyLowerTriangular);
190
 
      if((Mode == UpperTriangular && row <= col)
191
 
      || (Mode == LowerTriangular && row >= col)
192
 
      || (Mode == StrictlyUpperTriangular && row < col)
193
 
      || (Mode == StrictlyLowerTriangular && row > col))
194
 
        dst.copyCoeff(row, col, src);
195
 
    }
196
 
  }
197
 
};
198
 
 
199
 
template<typename Derived1, typename Derived2, unsigned int Mode>
200
 
struct ei_part_assignment_impl<Derived1, Derived2, Mode, 1>
201
 
{
202
 
  inline static void run(Derived1 &dst, const Derived2 &src)
203
 
  {
204
 
    if(!(Mode & ZeroDiagBit))
205
 
      dst.copyCoeff(0, 0, src);
206
 
  }
207
 
};
208
 
 
209
 
// prevent buggy user code from causing an infinite recursion
210
 
template<typename Derived1, typename Derived2, unsigned int Mode>
211
 
struct ei_part_assignment_impl<Derived1, Derived2, Mode, 0>
212
 
{
213
 
  inline static void run(Derived1 &, const Derived2 &) {}
214
 
};
215
 
 
216
 
template<typename Derived1, typename Derived2>
217
 
struct ei_part_assignment_impl<Derived1, Derived2, UpperTriangular, Dynamic>
218
 
{
219
 
  inline static void run(Derived1 &dst, const Derived2 &src)
220
 
  {
221
 
    for(int j = 0; j < dst.cols(); ++j)
222
 
      for(int i = 0; i <= j; ++i)
223
 
        dst.copyCoeff(i, j, src);
224
 
  }
225
 
};
226
 
 
227
 
template<typename Derived1, typename Derived2>
228
 
struct ei_part_assignment_impl<Derived1, Derived2, LowerTriangular, Dynamic>
229
 
{
230
 
  inline static void run(Derived1 &dst, const Derived2 &src)
231
 
  {
232
 
    for(int j = 0; j < dst.cols(); ++j)
233
 
      for(int i = j; i < dst.rows(); ++i)
234
 
        dst.copyCoeff(i, j, src);
235
 
  }
236
 
};
237
 
 
238
 
template<typename Derived1, typename Derived2>
239
 
struct ei_part_assignment_impl<Derived1, Derived2, StrictlyUpperTriangular, Dynamic>
240
 
{
241
 
  inline static void run(Derived1 &dst, const Derived2 &src)
242
 
  {
243
 
    for(int j = 0; j < dst.cols(); ++j)
244
 
      for(int i = 0; i < j; ++i)
245
 
        dst.copyCoeff(i, j, src);
246
 
  }
247
 
};
248
 
template<typename Derived1, typename Derived2>
249
 
struct ei_part_assignment_impl<Derived1, Derived2, StrictlyLowerTriangular, Dynamic>
250
 
{
251
 
  inline static void run(Derived1 &dst, const Derived2 &src)
252
 
  {
253
 
    for(int j = 0; j < dst.cols(); ++j)
254
 
      for(int i = j+1; i < dst.rows(); ++i)
255
 
        dst.copyCoeff(i, j, src);
256
 
  }
257
 
};
258
 
template<typename Derived1, typename Derived2>
259
 
struct ei_part_assignment_impl<Derived1, Derived2, SelfAdjoint, Dynamic>
260
 
{
261
 
  inline static void run(Derived1 &dst, const Derived2 &src)
262
 
  {
263
 
    for(int j = 0; j < dst.cols(); ++j)
264
 
    {
265
 
      for(int i = 0; i < j; ++i)
266
 
        dst.coeffRef(j, i) = ei_conj(dst.coeffRef(i, j) = src.coeff(i, j));
267
 
      dst.coeffRef(j, j) = ei_real(src.coeff(j, j));
268
 
    }
269
 
  }
270
 
};
271
 
 
272
 
template<typename MatrixType, unsigned int Mode>
273
 
template<typename Other>
274
 
void Part<MatrixType, Mode>::lazyAssign(const Other& other)
275
 
{
276
 
  const bool unroll = MatrixType::SizeAtCompileTime * Other::CoeffReadCost / 2 <= EIGEN_UNROLLING_LIMIT;
277
 
  ei_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols());
278
 
 
279
 
  ei_part_assignment_impl
280
 
    <MatrixType, Other, Mode,
281
 
    unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic
282
 
    >::run(m_matrix.const_cast_derived(), other.derived());
283
 
}
284
 
 
285
 
/** \nonstableyet
286
 
  * \returns a lvalue pseudo-expression allowing to perform special operations on \c *this.
287
 
  *
288
 
  * The \a Mode parameter can have the following values: \c UpperTriangular, \c StrictlyUpperTriangular, \c LowerTriangular,
289
 
  * \c StrictlyLowerTriangular, \c SelfAdjoint.
290
 
  *
291
 
  * \addexample PartExample \label How to write to a triangular part of a matrix
292
 
  *
293
 
  * Example: \include MatrixBase_part.cpp
294
 
  * Output: \verbinclude MatrixBase_part.out
295
 
  *
296
 
  * \sa class Part, MatrixBase::extract(), MatrixBase::marked()
297
 
  */
298
 
template<typename Derived>
299
 
template<unsigned int Mode>
300
 
inline Part<Derived, Mode> MatrixBase<Derived>::part()
301
 
{
302
 
  return Part<Derived, Mode>(derived());
303
 
}
304
 
 
305
 
/** \returns true if *this is approximately equal to an upper triangular matrix,
306
 
  *          within the precision given by \a prec.
307
 
  *
308
 
  * \sa isLowerTriangular(), extract(), part(), marked()
309
 
  */
310
 
template<typename Derived>
311
 
bool MatrixBase<Derived>::isUpperTriangular(RealScalar prec) const
312
 
{
313
 
  if(cols() != rows()) return false;
314
 
  RealScalar maxAbsOnUpperTriangularPart = static_cast<RealScalar>(-1);
315
 
  for(int j = 0; j < cols(); ++j)
316
 
    for(int i = 0; i <= j; ++i)
317
 
    {
318
 
      RealScalar absValue = ei_abs(coeff(i,j));
319
 
      if(absValue > maxAbsOnUpperTriangularPart) maxAbsOnUpperTriangularPart = absValue;
320
 
    }
321
 
  for(int j = 0; j < cols()-1; ++j)
322
 
    for(int i = j+1; i < rows(); ++i)
323
 
      if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnUpperTriangularPart, prec)) return false;
324
 
  return true;
325
 
}
326
 
 
327
 
/** \returns true if *this is approximately equal to a lower triangular matrix,
328
 
  *          within the precision given by \a prec.
329
 
  *
330
 
  * \sa isUpperTriangular(), extract(), part(), marked()
331
 
  */
332
 
template<typename Derived>
333
 
bool MatrixBase<Derived>::isLowerTriangular(RealScalar prec) const
334
 
{
335
 
  if(cols() != rows()) return false;
336
 
  RealScalar maxAbsOnLowerTriangularPart = static_cast<RealScalar>(-1);
337
 
  for(int j = 0; j < cols(); ++j)
338
 
    for(int i = j; i < rows(); ++i)
339
 
    {
340
 
      RealScalar absValue = ei_abs(coeff(i,j));
341
 
      if(absValue > maxAbsOnLowerTriangularPart) maxAbsOnLowerTriangularPart = absValue;
342
 
    }
343
 
  for(int j = 1; j < cols(); ++j)
344
 
    for(int i = 0; i < j; ++i)
345
 
      if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnLowerTriangularPart, prec)) return false;
346
 
  return true;
347
 
}
348
 
 
349
 
template<typename MatrixType, unsigned int Mode>
350
 
template<typename Other>
351
 
inline Part<MatrixType, Mode>& Part<MatrixType, Mode>::operator+=(const Other& other)
352
 
{
353
 
  return *this = m_matrix + other;
354
 
}
355
 
 
356
 
template<typename MatrixType, unsigned int Mode>
357
 
template<typename Other>
358
 
inline Part<MatrixType, Mode>& Part<MatrixType, Mode>::operator-=(const Other& other)
359
 
{
360
 
  return *this = m_matrix - other;
361
 
}
362
 
 
363
 
template<typename MatrixType, unsigned int Mode>
364
 
inline Part<MatrixType, Mode>& Part<MatrixType, Mode>::operator*=
365
 
(const typename ei_traits<MatrixType>::Scalar& other)
366
 
{
367
 
  return *this = m_matrix * other;
368
 
}
369
 
 
370
 
template<typename MatrixType, unsigned int Mode>
371
 
inline Part<MatrixType, Mode>& Part<MatrixType, Mode>::operator/=
372
 
(const typename ei_traits<MatrixType>::Scalar& other)
373
 
{
374
 
  return *this = m_matrix / other;
375
 
}
376
 
 
377
 
#endif // EIGEN_PART_H