1 | // This file is part of Eigen, a lightweight C++ template library
|
---|
2 | // for linear algebra.
|
---|
3 | //
|
---|
4 | // Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com>
|
---|
5 | // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
|
---|
6 | //
|
---|
7 | // This Source Code Form is subject to the terms of the Mozilla
|
---|
8 | // Public License v. 2.0. If a copy of the MPL was not distributed
|
---|
9 | // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
---|
10 |
|
---|
11 | #ifndef EIGEN_TRIANGULARMATRIX_H
|
---|
12 | #define EIGEN_TRIANGULARMATRIX_H
|
---|
13 |
|
---|
14 | namespace Eigen {
|
---|
15 |
|
---|
16 | namespace internal {
|
---|
17 |
|
---|
18 | template<int Side, typename TriangularType, typename Rhs> struct triangular_solve_retval;
|
---|
19 |
|
---|
20 | }
|
---|
21 |
|
---|
22 | /** \internal
|
---|
23 | *
|
---|
24 | * \class TriangularBase
|
---|
25 | * \ingroup Core_Module
|
---|
26 | *
|
---|
27 | * \brief Base class for triangular part in a matrix
|
---|
28 | */
|
---|
29 | template<typename Derived> class TriangularBase : public EigenBase<Derived>
|
---|
30 | {
|
---|
31 | public:
|
---|
32 |
|
---|
33 | enum {
|
---|
34 | Mode = internal::traits<Derived>::Mode,
|
---|
35 | CoeffReadCost = internal::traits<Derived>::CoeffReadCost,
|
---|
36 | RowsAtCompileTime = internal::traits<Derived>::RowsAtCompileTime,
|
---|
37 | ColsAtCompileTime = internal::traits<Derived>::ColsAtCompileTime,
|
---|
38 | MaxRowsAtCompileTime = internal::traits<Derived>::MaxRowsAtCompileTime,
|
---|
39 | MaxColsAtCompileTime = internal::traits<Derived>::MaxColsAtCompileTime
|
---|
40 | };
|
---|
41 | typedef typename internal::traits<Derived>::Scalar Scalar;
|
---|
42 | typedef typename internal::traits<Derived>::StorageKind StorageKind;
|
---|
43 | typedef typename internal::traits<Derived>::Index Index;
|
---|
44 | typedef typename internal::traits<Derived>::DenseMatrixType DenseMatrixType;
|
---|
45 | typedef DenseMatrixType DenseType;
|
---|
46 |
|
---|
47 | inline TriangularBase() { eigen_assert(!((Mode&UnitDiag) && (Mode&ZeroDiag))); }
|
---|
48 |
|
---|
49 | inline Index rows() const { return derived().rows(); }
|
---|
50 | inline Index cols() const { return derived().cols(); }
|
---|
51 | inline Index outerStride() const { return derived().outerStride(); }
|
---|
52 | inline Index innerStride() const { return derived().innerStride(); }
|
---|
53 |
|
---|
54 | inline Scalar coeff(Index row, Index col) const { return derived().coeff(row,col); }
|
---|
55 | inline Scalar& coeffRef(Index row, Index col) { return derived().coeffRef(row,col); }
|
---|
56 |
|
---|
57 | /** \see MatrixBase::copyCoeff(row,col)
|
---|
58 | */
|
---|
59 | template<typename Other>
|
---|
60 | EIGEN_STRONG_INLINE void copyCoeff(Index row, Index col, Other& other)
|
---|
61 | {
|
---|
62 | derived().coeffRef(row, col) = other.coeff(row, col);
|
---|
63 | }
|
---|
64 |
|
---|
65 | inline Scalar operator()(Index row, Index col) const
|
---|
66 | {
|
---|
67 | check_coordinates(row, col);
|
---|
68 | return coeff(row,col);
|
---|
69 | }
|
---|
70 | inline Scalar& operator()(Index row, Index col)
|
---|
71 | {
|
---|
72 | check_coordinates(row, col);
|
---|
73 | return coeffRef(row,col);
|
---|
74 | }
|
---|
75 |
|
---|
76 | #ifndef EIGEN_PARSED_BY_DOXYGEN
|
---|
77 | inline const Derived& derived() const { return *static_cast<const Derived*>(this); }
|
---|
78 | inline Derived& derived() { return *static_cast<Derived*>(this); }
|
---|
79 | #endif // not EIGEN_PARSED_BY_DOXYGEN
|
---|
80 |
|
---|
81 | template<typename DenseDerived>
|
---|
82 | void evalTo(MatrixBase<DenseDerived> &other) const;
|
---|
83 | template<typename DenseDerived>
|
---|
84 | void evalToLazy(MatrixBase<DenseDerived> &other) const;
|
---|
85 |
|
---|
86 | DenseMatrixType toDenseMatrix() const
|
---|
87 | {
|
---|
88 | DenseMatrixType res(rows(), cols());
|
---|
89 | evalToLazy(res);
|
---|
90 | return res;
|
---|
91 | }
|
---|
92 |
|
---|
93 | protected:
|
---|
94 |
|
---|
95 | void check_coordinates(Index row, Index col) const
|
---|
96 | {
|
---|
97 | EIGEN_ONLY_USED_FOR_DEBUG(row);
|
---|
98 | EIGEN_ONLY_USED_FOR_DEBUG(col);
|
---|
99 | eigen_assert(col>=0 && col<cols() && row>=0 && row<rows());
|
---|
100 | const int mode = int(Mode) & ~SelfAdjoint;
|
---|
101 | EIGEN_ONLY_USED_FOR_DEBUG(mode);
|
---|
102 | eigen_assert((mode==Upper && col>=row)
|
---|
103 | || (mode==Lower && col<=row)
|
---|
104 | || ((mode==StrictlyUpper || mode==UnitUpper) && col>row)
|
---|
105 | || ((mode==StrictlyLower || mode==UnitLower) && col<row));
|
---|
106 | }
|
---|
107 |
|
---|
108 | #ifdef EIGEN_INTERNAL_DEBUGGING
|
---|
109 | void check_coordinates_internal(Index row, Index col) const
|
---|
110 | {
|
---|
111 | check_coordinates(row, col);
|
---|
112 | }
|
---|
113 | #else
|
---|
114 | void check_coordinates_internal(Index , Index ) const {}
|
---|
115 | #endif
|
---|
116 |
|
---|
117 | };
|
---|
118 |
|
---|
119 | /** \class TriangularView
|
---|
120 | * \ingroup Core_Module
|
---|
121 | *
|
---|
122 | * \brief Base class for triangular part in a matrix
|
---|
123 | *
|
---|
124 | * \param MatrixType the type of the object in which we are taking the triangular part
|
---|
125 | * \param Mode the kind of triangular matrix expression to construct. Can be #Upper,
|
---|
126 | * #Lower, #UnitUpper, #UnitLower, #StrictlyUpper, or #StrictlyLower.
|
---|
127 | * This is in fact a bit field; it must have either #Upper or #Lower,
|
---|
128 | * and additionnaly it may have #UnitDiag or #ZeroDiag or neither.
|
---|
129 | *
|
---|
130 | * This class represents a triangular part of a matrix, not necessarily square. Strictly speaking, for rectangular
|
---|
131 | * matrices one should speak of "trapezoid" parts. This class is the return type
|
---|
132 | * of MatrixBase::triangularView() and most of the time this is the only way it is used.
|
---|
133 | *
|
---|
134 | * \sa MatrixBase::triangularView()
|
---|
135 | */
|
---|
136 | namespace internal {
|
---|
137 | template<typename MatrixType, unsigned int _Mode>
|
---|
138 | struct traits<TriangularView<MatrixType, _Mode> > : traits<MatrixType>
|
---|
139 | {
|
---|
140 | typedef typename nested<MatrixType>::type MatrixTypeNested;
|
---|
141 | typedef typename remove_reference<MatrixTypeNested>::type MatrixTypeNestedNonRef;
|
---|
142 | typedef typename remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned;
|
---|
143 | typedef MatrixType ExpressionType;
|
---|
144 | typedef typename MatrixType::PlainObject DenseMatrixType;
|
---|
145 | enum {
|
---|
146 | Mode = _Mode,
|
---|
147 | Flags = (MatrixTypeNestedCleaned::Flags & (HereditaryBits) & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit))) | Mode,
|
---|
148 | CoeffReadCost = MatrixTypeNestedCleaned::CoeffReadCost
|
---|
149 | };
|
---|
150 | };
|
---|
151 | }
|
---|
152 |
|
---|
153 | template<int Mode, bool LhsIsTriangular,
|
---|
154 | typename Lhs, bool LhsIsVector,
|
---|
155 | typename Rhs, bool RhsIsVector>
|
---|
156 | struct TriangularProduct;
|
---|
157 |
|
---|
158 | template<typename _MatrixType, unsigned int _Mode> class TriangularView
|
---|
159 | : public TriangularBase<TriangularView<_MatrixType, _Mode> >
|
---|
160 | {
|
---|
161 | public:
|
---|
162 |
|
---|
163 | typedef TriangularBase<TriangularView> Base;
|
---|
164 | typedef typename internal::traits<TriangularView>::Scalar Scalar;
|
---|
165 |
|
---|
166 | typedef _MatrixType MatrixType;
|
---|
167 | typedef typename internal::traits<TriangularView>::DenseMatrixType DenseMatrixType;
|
---|
168 | typedef DenseMatrixType PlainObject;
|
---|
169 |
|
---|
170 | protected:
|
---|
171 | typedef typename internal::traits<TriangularView>::MatrixTypeNested MatrixTypeNested;
|
---|
172 | typedef typename internal::traits<TriangularView>::MatrixTypeNestedNonRef MatrixTypeNestedNonRef;
|
---|
173 | typedef typename internal::traits<TriangularView>::MatrixTypeNestedCleaned MatrixTypeNestedCleaned;
|
---|
174 |
|
---|
175 | typedef typename internal::remove_all<typename MatrixType::ConjugateReturnType>::type MatrixConjugateReturnType;
|
---|
176 |
|
---|
177 | public:
|
---|
178 | using Base::evalToLazy;
|
---|
179 |
|
---|
180 |
|
---|
181 | typedef typename internal::traits<TriangularView>::StorageKind StorageKind;
|
---|
182 | typedef typename internal::traits<TriangularView>::Index Index;
|
---|
183 |
|
---|
184 | enum {
|
---|
185 | Mode = _Mode,
|
---|
186 | TransposeMode = (Mode & Upper ? Lower : 0)
|
---|
187 | | (Mode & Lower ? Upper : 0)
|
---|
188 | | (Mode & (UnitDiag))
|
---|
189 | | (Mode & (ZeroDiag))
|
---|
190 | };
|
---|
191 |
|
---|
192 | inline TriangularView(const MatrixType& matrix) : m_matrix(matrix)
|
---|
193 | {}
|
---|
194 |
|
---|
195 | inline Index rows() const { return m_matrix.rows(); }
|
---|
196 | inline Index cols() const { return m_matrix.cols(); }
|
---|
197 | inline Index outerStride() const { return m_matrix.outerStride(); }
|
---|
198 | inline Index innerStride() const { return m_matrix.innerStride(); }
|
---|
199 |
|
---|
200 | /** \sa MatrixBase::operator+=() */
|
---|
201 | template<typename Other> TriangularView& operator+=(const DenseBase<Other>& other) { return *this = m_matrix + other.derived(); }
|
---|
202 | /** \sa MatrixBase::operator-=() */
|
---|
203 | template<typename Other> TriangularView& operator-=(const DenseBase<Other>& other) { return *this = m_matrix - other.derived(); }
|
---|
204 | /** \sa MatrixBase::operator*=() */
|
---|
205 | TriangularView& operator*=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = m_matrix * other; }
|
---|
206 | /** \sa MatrixBase::operator/=() */
|
---|
207 | TriangularView& operator/=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = m_matrix / other; }
|
---|
208 |
|
---|
209 | /** \sa MatrixBase::fill() */
|
---|
210 | void fill(const Scalar& value) { setConstant(value); }
|
---|
211 | /** \sa MatrixBase::setConstant() */
|
---|
212 | TriangularView& setConstant(const Scalar& value)
|
---|
213 | { return *this = MatrixType::Constant(rows(), cols(), value); }
|
---|
214 | /** \sa MatrixBase::setZero() */
|
---|
215 | TriangularView& setZero() { return setConstant(Scalar(0)); }
|
---|
216 | /** \sa MatrixBase::setOnes() */
|
---|
217 | TriangularView& setOnes() { return setConstant(Scalar(1)); }
|
---|
218 |
|
---|
219 | /** \sa MatrixBase::coeff()
|
---|
220 | * \warning the coordinates must fit into the referenced triangular part
|
---|
221 | */
|
---|
222 | inline Scalar coeff(Index row, Index col) const
|
---|
223 | {
|
---|
224 | Base::check_coordinates_internal(row, col);
|
---|
225 | return m_matrix.coeff(row, col);
|
---|
226 | }
|
---|
227 |
|
---|
228 | /** \sa MatrixBase::coeffRef()
|
---|
229 | * \warning the coordinates must fit into the referenced triangular part
|
---|
230 | */
|
---|
231 | inline Scalar& coeffRef(Index row, Index col)
|
---|
232 | {
|
---|
233 | Base::check_coordinates_internal(row, col);
|
---|
234 | return m_matrix.const_cast_derived().coeffRef(row, col);
|
---|
235 | }
|
---|
236 |
|
---|
237 | const MatrixTypeNestedCleaned& nestedExpression() const { return m_matrix; }
|
---|
238 | MatrixTypeNestedCleaned& nestedExpression() { return *const_cast<MatrixTypeNestedCleaned*>(&m_matrix); }
|
---|
239 |
|
---|
240 | /** Assigns a triangular matrix to a triangular part of a dense matrix */
|
---|
241 | template<typename OtherDerived>
|
---|
242 | TriangularView& operator=(const TriangularBase<OtherDerived>& other);
|
---|
243 |
|
---|
244 | template<typename OtherDerived>
|
---|
245 | TriangularView& operator=(const MatrixBase<OtherDerived>& other);
|
---|
246 |
|
---|
247 | TriangularView& operator=(const TriangularView& other)
|
---|
248 | { return *this = other.nestedExpression(); }
|
---|
249 |
|
---|
250 | template<typename OtherDerived>
|
---|
251 | void lazyAssign(const TriangularBase<OtherDerived>& other);
|
---|
252 |
|
---|
253 | template<typename OtherDerived>
|
---|
254 | void lazyAssign(const MatrixBase<OtherDerived>& other);
|
---|
255 |
|
---|
256 | /** \sa MatrixBase::conjugate() */
|
---|
257 | inline TriangularView<MatrixConjugateReturnType,Mode> conjugate()
|
---|
258 | { return m_matrix.conjugate(); }
|
---|
259 | /** \sa MatrixBase::conjugate() const */
|
---|
260 | inline const TriangularView<MatrixConjugateReturnType,Mode> conjugate() const
|
---|
261 | { return m_matrix.conjugate(); }
|
---|
262 |
|
---|
263 | /** \sa MatrixBase::adjoint() const */
|
---|
264 | inline const TriangularView<const typename MatrixType::AdjointReturnType,TransposeMode> adjoint() const
|
---|
265 | { return m_matrix.adjoint(); }
|
---|
266 |
|
---|
267 | /** \sa MatrixBase::transpose() */
|
---|
268 | inline TriangularView<Transpose<MatrixType>,TransposeMode> transpose()
|
---|
269 | {
|
---|
270 | EIGEN_STATIC_ASSERT_LVALUE(MatrixType)
|
---|
271 | return m_matrix.const_cast_derived().transpose();
|
---|
272 | }
|
---|
273 | /** \sa MatrixBase::transpose() const */
|
---|
274 | inline const TriangularView<Transpose<MatrixType>,TransposeMode> transpose() const
|
---|
275 | {
|
---|
276 | return m_matrix.transpose();
|
---|
277 | }
|
---|
278 |
|
---|
279 | /** Efficient triangular matrix times vector/matrix product */
|
---|
280 | template<typename OtherDerived>
|
---|
281 | TriangularProduct<Mode, true, MatrixType, false, OtherDerived, OtherDerived::ColsAtCompileTime==1>
|
---|
282 | operator*(const MatrixBase<OtherDerived>& rhs) const
|
---|
283 | {
|
---|
284 | return TriangularProduct
|
---|
285 | <Mode, true, MatrixType, false, OtherDerived, OtherDerived::ColsAtCompileTime==1>
|
---|
286 | (m_matrix, rhs.derived());
|
---|
287 | }
|
---|
288 |
|
---|
289 | /** Efficient vector/matrix times triangular matrix product */
|
---|
290 | template<typename OtherDerived> friend
|
---|
291 | TriangularProduct<Mode, false, OtherDerived, OtherDerived::RowsAtCompileTime==1, MatrixType, false>
|
---|
292 | operator*(const MatrixBase<OtherDerived>& lhs, const TriangularView& rhs)
|
---|
293 | {
|
---|
294 | return TriangularProduct
|
---|
295 | <Mode, false, OtherDerived, OtherDerived::RowsAtCompileTime==1, MatrixType, false>
|
---|
296 | (lhs.derived(),rhs.m_matrix);
|
---|
297 | }
|
---|
298 |
|
---|
299 | #ifdef EIGEN2_SUPPORT
|
---|
300 | template<typename OtherDerived>
|
---|
301 | struct eigen2_product_return_type
|
---|
302 | {
|
---|
303 | typedef typename TriangularView<MatrixType,Mode>::DenseMatrixType DenseMatrixType;
|
---|
304 | typedef typename OtherDerived::PlainObject::DenseType OtherPlainObject;
|
---|
305 | typedef typename ProductReturnType<DenseMatrixType, OtherPlainObject>::Type ProdRetType;
|
---|
306 | typedef typename ProdRetType::PlainObject type;
|
---|
307 | };
|
---|
308 | template<typename OtherDerived>
|
---|
309 | const typename eigen2_product_return_type<OtherDerived>::type
|
---|
310 | operator*(const EigenBase<OtherDerived>& rhs) const
|
---|
311 | {
|
---|
312 | typename OtherDerived::PlainObject::DenseType rhsPlainObject;
|
---|
313 | rhs.evalTo(rhsPlainObject);
|
---|
314 | return this->toDenseMatrix() * rhsPlainObject;
|
---|
315 | }
|
---|
316 | template<typename OtherMatrixType>
|
---|
317 | bool isApprox(const TriangularView<OtherMatrixType, Mode>& other, typename NumTraits<Scalar>::Real precision = NumTraits<Scalar>::dummy_precision()) const
|
---|
318 | {
|
---|
319 | return this->toDenseMatrix().isApprox(other.toDenseMatrix(), precision);
|
---|
320 | }
|
---|
321 | template<typename OtherDerived>
|
---|
322 | bool isApprox(const MatrixBase<OtherDerived>& other, typename NumTraits<Scalar>::Real precision = NumTraits<Scalar>::dummy_precision()) const
|
---|
323 | {
|
---|
324 | return this->toDenseMatrix().isApprox(other, precision);
|
---|
325 | }
|
---|
326 | #endif // EIGEN2_SUPPORT
|
---|
327 |
|
---|
328 | template<int Side, typename Other>
|
---|
329 | inline const internal::triangular_solve_retval<Side,TriangularView, Other>
|
---|
330 | solve(const MatrixBase<Other>& other) const;
|
---|
331 |
|
---|
332 | template<int Side, typename OtherDerived>
|
---|
333 | void solveInPlace(const MatrixBase<OtherDerived>& other) const;
|
---|
334 |
|
---|
335 | template<typename Other>
|
---|
336 | inline const internal::triangular_solve_retval<OnTheLeft,TriangularView, Other>
|
---|
337 | solve(const MatrixBase<Other>& other) const
|
---|
338 | { return solve<OnTheLeft>(other); }
|
---|
339 |
|
---|
340 | template<typename OtherDerived>
|
---|
341 | void solveInPlace(const MatrixBase<OtherDerived>& other) const
|
---|
342 | { return solveInPlace<OnTheLeft>(other); }
|
---|
343 |
|
---|
344 | const SelfAdjointView<MatrixTypeNestedNonRef,Mode> selfadjointView() const
|
---|
345 | {
|
---|
346 | EIGEN_STATIC_ASSERT((Mode&UnitDiag)==0,PROGRAMMING_ERROR);
|
---|
347 | return SelfAdjointView<MatrixTypeNestedNonRef,Mode>(m_matrix);
|
---|
348 | }
|
---|
349 | SelfAdjointView<MatrixTypeNestedNonRef,Mode> selfadjointView()
|
---|
350 | {
|
---|
351 | EIGEN_STATIC_ASSERT((Mode&UnitDiag)==0,PROGRAMMING_ERROR);
|
---|
352 | return SelfAdjointView<MatrixTypeNestedNonRef,Mode>(m_matrix);
|
---|
353 | }
|
---|
354 |
|
---|
355 | template<typename OtherDerived>
|
---|
356 | void swap(TriangularBase<OtherDerived> const & other)
|
---|
357 | {
|
---|
358 | TriangularView<SwapWrapper<MatrixType>,Mode>(const_cast<MatrixType&>(m_matrix)).lazyAssign(other.derived());
|
---|
359 | }
|
---|
360 |
|
---|
361 | template<typename OtherDerived>
|
---|
362 | void swap(MatrixBase<OtherDerived> const & other)
|
---|
363 | {
|
---|
364 | SwapWrapper<MatrixType> swaper(const_cast<MatrixType&>(m_matrix));
|
---|
365 | TriangularView<SwapWrapper<MatrixType>,Mode>(swaper).lazyAssign(other.derived());
|
---|
366 | }
|
---|
367 |
|
---|
368 | Scalar determinant() const
|
---|
369 | {
|
---|
370 | if (Mode & UnitDiag)
|
---|
371 | return 1;
|
---|
372 | else if (Mode & ZeroDiag)
|
---|
373 | return 0;
|
---|
374 | else
|
---|
375 | return m_matrix.diagonal().prod();
|
---|
376 | }
|
---|
377 |
|
---|
378 | // TODO simplify the following:
|
---|
379 | template<typename ProductDerived, typename Lhs, typename Rhs>
|
---|
380 | EIGEN_STRONG_INLINE TriangularView& operator=(const ProductBase<ProductDerived, Lhs,Rhs>& other)
|
---|
381 | {
|
---|
382 | setZero();
|
---|
383 | return assignProduct(other.derived(),1);
|
---|
384 | }
|
---|
385 |
|
---|
386 | template<typename ProductDerived, typename Lhs, typename Rhs>
|
---|
387 | EIGEN_STRONG_INLINE TriangularView& operator+=(const ProductBase<ProductDerived, Lhs,Rhs>& other)
|
---|
388 | {
|
---|
389 | return assignProduct(other.derived(),1);
|
---|
390 | }
|
---|
391 |
|
---|
392 | template<typename ProductDerived, typename Lhs, typename Rhs>
|
---|
393 | EIGEN_STRONG_INLINE TriangularView& operator-=(const ProductBase<ProductDerived, Lhs,Rhs>& other)
|
---|
394 | {
|
---|
395 | return assignProduct(other.derived(),-1);
|
---|
396 | }
|
---|
397 |
|
---|
398 |
|
---|
399 | template<typename ProductDerived>
|
---|
400 | EIGEN_STRONG_INLINE TriangularView& operator=(const ScaledProduct<ProductDerived>& other)
|
---|
401 | {
|
---|
402 | setZero();
|
---|
403 | return assignProduct(other.derived(),other.alpha());
|
---|
404 | }
|
---|
405 |
|
---|
406 | template<typename ProductDerived>
|
---|
407 | EIGEN_STRONG_INLINE TriangularView& operator+=(const ScaledProduct<ProductDerived>& other)
|
---|
408 | {
|
---|
409 | return assignProduct(other.derived(),other.alpha());
|
---|
410 | }
|
---|
411 |
|
---|
412 | template<typename ProductDerived>
|
---|
413 | EIGEN_STRONG_INLINE TriangularView& operator-=(const ScaledProduct<ProductDerived>& other)
|
---|
414 | {
|
---|
415 | return assignProduct(other.derived(),-other.alpha());
|
---|
416 | }
|
---|
417 |
|
---|
418 | protected:
|
---|
419 |
|
---|
420 | template<typename ProductDerived, typename Lhs, typename Rhs>
|
---|
421 | EIGEN_STRONG_INLINE TriangularView& assignProduct(const ProductBase<ProductDerived, Lhs,Rhs>& prod, const Scalar& alpha);
|
---|
422 |
|
---|
423 | template<int Mode, bool LhsIsTriangular,
|
---|
424 | typename Lhs, bool LhsIsVector,
|
---|
425 | typename Rhs, bool RhsIsVector>
|
---|
426 | EIGEN_STRONG_INLINE TriangularView& assignProduct(const TriangularProduct<Mode, LhsIsTriangular, Lhs, LhsIsVector, Rhs, RhsIsVector>& prod, const Scalar& alpha)
|
---|
427 | {
|
---|
428 | lazyAssign(alpha*prod.eval());
|
---|
429 | return *this;
|
---|
430 | }
|
---|
431 |
|
---|
432 | MatrixTypeNested m_matrix;
|
---|
433 | };
|
---|
434 |
|
---|
435 | /***************************************************************************
|
---|
436 | * Implementation of triangular evaluation/assignment
|
---|
437 | ***************************************************************************/
|
---|
438 |
|
---|
439 | namespace internal {
|
---|
440 |
|
---|
441 | template<typename Derived1, typename Derived2, unsigned int Mode, int UnrollCount, bool ClearOpposite>
|
---|
442 | struct triangular_assignment_selector
|
---|
443 | {
|
---|
444 | enum {
|
---|
445 | col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
|
---|
446 | row = (UnrollCount-1) % Derived1::RowsAtCompileTime
|
---|
447 | };
|
---|
448 |
|
---|
449 | typedef typename Derived1::Scalar Scalar;
|
---|
450 |
|
---|
451 | static inline void run(Derived1 &dst, const Derived2 &src)
|
---|
452 | {
|
---|
453 | triangular_assignment_selector<Derived1, Derived2, Mode, UnrollCount-1, ClearOpposite>::run(dst, src);
|
---|
454 |
|
---|
455 | eigen_assert( Mode == Upper || Mode == Lower
|
---|
456 | || Mode == StrictlyUpper || Mode == StrictlyLower
|
---|
457 | || Mode == UnitUpper || Mode == UnitLower);
|
---|
458 | if((Mode == Upper && row <= col)
|
---|
459 | || (Mode == Lower && row >= col)
|
---|
460 | || (Mode == StrictlyUpper && row < col)
|
---|
461 | || (Mode == StrictlyLower && row > col)
|
---|
462 | || (Mode == UnitUpper && row < col)
|
---|
463 | || (Mode == UnitLower && row > col))
|
---|
464 | dst.copyCoeff(row, col, src);
|
---|
465 | else if(ClearOpposite)
|
---|
466 | {
|
---|
467 | if (Mode&UnitDiag && row==col)
|
---|
468 | dst.coeffRef(row, col) = Scalar(1);
|
---|
469 | else
|
---|
470 | dst.coeffRef(row, col) = Scalar(0);
|
---|
471 | }
|
---|
472 | }
|
---|
473 | };
|
---|
474 |
|
---|
475 | // prevent buggy user code from causing an infinite recursion
|
---|
476 | template<typename Derived1, typename Derived2, unsigned int Mode, bool ClearOpposite>
|
---|
477 | struct triangular_assignment_selector<Derived1, Derived2, Mode, 0, ClearOpposite>
|
---|
478 | {
|
---|
479 | static inline void run(Derived1 &, const Derived2 &) {}
|
---|
480 | };
|
---|
481 |
|
---|
482 | template<typename Derived1, typename Derived2, bool ClearOpposite>
|
---|
483 | struct triangular_assignment_selector<Derived1, Derived2, Upper, Dynamic, ClearOpposite>
|
---|
484 | {
|
---|
485 | typedef typename Derived1::Index Index;
|
---|
486 | typedef typename Derived1::Scalar Scalar;
|
---|
487 | static inline void run(Derived1 &dst, const Derived2 &src)
|
---|
488 | {
|
---|
489 | for(Index j = 0; j < dst.cols(); ++j)
|
---|
490 | {
|
---|
491 | Index maxi = (std::min)(j, dst.rows()-1);
|
---|
492 | for(Index i = 0; i <= maxi; ++i)
|
---|
493 | dst.copyCoeff(i, j, src);
|
---|
494 | if (ClearOpposite)
|
---|
495 | for(Index i = maxi+1; i < dst.rows(); ++i)
|
---|
496 | dst.coeffRef(i, j) = Scalar(0);
|
---|
497 | }
|
---|
498 | }
|
---|
499 | };
|
---|
500 |
|
---|
501 | template<typename Derived1, typename Derived2, bool ClearOpposite>
|
---|
502 | struct triangular_assignment_selector<Derived1, Derived2, Lower, Dynamic, ClearOpposite>
|
---|
503 | {
|
---|
504 | typedef typename Derived1::Index Index;
|
---|
505 | static inline void run(Derived1 &dst, const Derived2 &src)
|
---|
506 | {
|
---|
507 | for(Index j = 0; j < dst.cols(); ++j)
|
---|
508 | {
|
---|
509 | for(Index i = j; i < dst.rows(); ++i)
|
---|
510 | dst.copyCoeff(i, j, src);
|
---|
511 | Index maxi = (std::min)(j, dst.rows());
|
---|
512 | if (ClearOpposite)
|
---|
513 | for(Index i = 0; i < maxi; ++i)
|
---|
514 | dst.coeffRef(i, j) = static_cast<typename Derived1::Scalar>(0);
|
---|
515 | }
|
---|
516 | }
|
---|
517 | };
|
---|
518 |
|
---|
519 | template<typename Derived1, typename Derived2, bool ClearOpposite>
|
---|
520 | struct triangular_assignment_selector<Derived1, Derived2, StrictlyUpper, Dynamic, ClearOpposite>
|
---|
521 | {
|
---|
522 | typedef typename Derived1::Index Index;
|
---|
523 | typedef typename Derived1::Scalar Scalar;
|
---|
524 | static inline void run(Derived1 &dst, const Derived2 &src)
|
---|
525 | {
|
---|
526 | for(Index j = 0; j < dst.cols(); ++j)
|
---|
527 | {
|
---|
528 | Index maxi = (std::min)(j, dst.rows());
|
---|
529 | for(Index i = 0; i < maxi; ++i)
|
---|
530 | dst.copyCoeff(i, j, src);
|
---|
531 | if (ClearOpposite)
|
---|
532 | for(Index i = maxi; i < dst.rows(); ++i)
|
---|
533 | dst.coeffRef(i, j) = Scalar(0);
|
---|
534 | }
|
---|
535 | }
|
---|
536 | };
|
---|
537 |
|
---|
538 | template<typename Derived1, typename Derived2, bool ClearOpposite>
|
---|
539 | struct triangular_assignment_selector<Derived1, Derived2, StrictlyLower, Dynamic, ClearOpposite>
|
---|
540 | {
|
---|
541 | typedef typename Derived1::Index Index;
|
---|
542 | static inline void run(Derived1 &dst, const Derived2 &src)
|
---|
543 | {
|
---|
544 | for(Index j = 0; j < dst.cols(); ++j)
|
---|
545 | {
|
---|
546 | for(Index i = j+1; i < dst.rows(); ++i)
|
---|
547 | dst.copyCoeff(i, j, src);
|
---|
548 | Index maxi = (std::min)(j, dst.rows()-1);
|
---|
549 | if (ClearOpposite)
|
---|
550 | for(Index i = 0; i <= maxi; ++i)
|
---|
551 | dst.coeffRef(i, j) = static_cast<typename Derived1::Scalar>(0);
|
---|
552 | }
|
---|
553 | }
|
---|
554 | };
|
---|
555 |
|
---|
556 | template<typename Derived1, typename Derived2, bool ClearOpposite>
|
---|
557 | struct triangular_assignment_selector<Derived1, Derived2, UnitUpper, Dynamic, ClearOpposite>
|
---|
558 | {
|
---|
559 | typedef typename Derived1::Index Index;
|
---|
560 | static inline void run(Derived1 &dst, const Derived2 &src)
|
---|
561 | {
|
---|
562 | for(Index j = 0; j < dst.cols(); ++j)
|
---|
563 | {
|
---|
564 | Index maxi = (std::min)(j, dst.rows());
|
---|
565 | for(Index i = 0; i < maxi; ++i)
|
---|
566 | dst.copyCoeff(i, j, src);
|
---|
567 | if (ClearOpposite)
|
---|
568 | {
|
---|
569 | for(Index i = maxi+1; i < dst.rows(); ++i)
|
---|
570 | dst.coeffRef(i, j) = 0;
|
---|
571 | }
|
---|
572 | }
|
---|
573 | dst.diagonal().setOnes();
|
---|
574 | }
|
---|
575 | };
|
---|
576 | template<typename Derived1, typename Derived2, bool ClearOpposite>
|
---|
577 | struct triangular_assignment_selector<Derived1, Derived2, UnitLower, Dynamic, ClearOpposite>
|
---|
578 | {
|
---|
579 | typedef typename Derived1::Index Index;
|
---|
580 | static inline void run(Derived1 &dst, const Derived2 &src)
|
---|
581 | {
|
---|
582 | for(Index j = 0; j < dst.cols(); ++j)
|
---|
583 | {
|
---|
584 | Index maxi = (std::min)(j, dst.rows());
|
---|
585 | for(Index i = maxi+1; i < dst.rows(); ++i)
|
---|
586 | dst.copyCoeff(i, j, src);
|
---|
587 | if (ClearOpposite)
|
---|
588 | {
|
---|
589 | for(Index i = 0; i < maxi; ++i)
|
---|
590 | dst.coeffRef(i, j) = 0;
|
---|
591 | }
|
---|
592 | }
|
---|
593 | dst.diagonal().setOnes();
|
---|
594 | }
|
---|
595 | };
|
---|
596 |
|
---|
597 | } // end namespace internal
|
---|
598 |
|
---|
599 | // FIXME should we keep that possibility
|
---|
600 | template<typename MatrixType, unsigned int Mode>
|
---|
601 | template<typename OtherDerived>
|
---|
602 | inline TriangularView<MatrixType, Mode>&
|
---|
603 | TriangularView<MatrixType, Mode>::operator=(const MatrixBase<OtherDerived>& other)
|
---|
604 | {
|
---|
605 | if(OtherDerived::Flags & EvalBeforeAssigningBit)
|
---|
606 | {
|
---|
607 | typename internal::plain_matrix_type<OtherDerived>::type other_evaluated(other.rows(), other.cols());
|
---|
608 | other_evaluated.template triangularView<Mode>().lazyAssign(other.derived());
|
---|
609 | lazyAssign(other_evaluated);
|
---|
610 | }
|
---|
611 | else
|
---|
612 | lazyAssign(other.derived());
|
---|
613 | return *this;
|
---|
614 | }
|
---|
615 |
|
---|
616 | // FIXME should we keep that possibility
|
---|
617 | template<typename MatrixType, unsigned int Mode>
|
---|
618 | template<typename OtherDerived>
|
---|
619 | void TriangularView<MatrixType, Mode>::lazyAssign(const MatrixBase<OtherDerived>& other)
|
---|
620 | {
|
---|
621 | enum {
|
---|
622 | unroll = MatrixType::SizeAtCompileTime != Dynamic
|
---|
623 | && internal::traits<OtherDerived>::CoeffReadCost != Dynamic
|
---|
624 | && MatrixType::SizeAtCompileTime*internal::traits<OtherDerived>::CoeffReadCost/2 <= EIGEN_UNROLLING_LIMIT
|
---|
625 | };
|
---|
626 | eigen_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols());
|
---|
627 |
|
---|
628 | internal::triangular_assignment_selector
|
---|
629 | <MatrixType, OtherDerived, int(Mode),
|
---|
630 | unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic,
|
---|
631 | false // do not change the opposite triangular part
|
---|
632 | >::run(m_matrix.const_cast_derived(), other.derived());
|
---|
633 | }
|
---|
634 |
|
---|
635 |
|
---|
636 |
|
---|
637 | template<typename MatrixType, unsigned int Mode>
|
---|
638 | template<typename OtherDerived>
|
---|
639 | inline TriangularView<MatrixType, Mode>&
|
---|
640 | TriangularView<MatrixType, Mode>::operator=(const TriangularBase<OtherDerived>& other)
|
---|
641 | {
|
---|
642 | eigen_assert(Mode == int(OtherDerived::Mode));
|
---|
643 | if(internal::traits<OtherDerived>::Flags & EvalBeforeAssigningBit)
|
---|
644 | {
|
---|
645 | typename OtherDerived::DenseMatrixType other_evaluated(other.rows(), other.cols());
|
---|
646 | other_evaluated.template triangularView<Mode>().lazyAssign(other.derived().nestedExpression());
|
---|
647 | lazyAssign(other_evaluated);
|
---|
648 | }
|
---|
649 | else
|
---|
650 | lazyAssign(other.derived().nestedExpression());
|
---|
651 | return *this;
|
---|
652 | }
|
---|
653 |
|
---|
654 | template<typename MatrixType, unsigned int Mode>
|
---|
655 | template<typename OtherDerived>
|
---|
656 | void TriangularView<MatrixType, Mode>::lazyAssign(const TriangularBase<OtherDerived>& other)
|
---|
657 | {
|
---|
658 | enum {
|
---|
659 | unroll = MatrixType::SizeAtCompileTime != Dynamic
|
---|
660 | && internal::traits<OtherDerived>::CoeffReadCost != Dynamic
|
---|
661 | && MatrixType::SizeAtCompileTime * internal::traits<OtherDerived>::CoeffReadCost / 2
|
---|
662 | <= EIGEN_UNROLLING_LIMIT
|
---|
663 | };
|
---|
664 | eigen_assert(m_matrix.rows() == other.rows() && m_matrix.cols() == other.cols());
|
---|
665 |
|
---|
666 | internal::triangular_assignment_selector
|
---|
667 | <MatrixType, OtherDerived, int(Mode),
|
---|
668 | unroll ? int(MatrixType::SizeAtCompileTime) : Dynamic,
|
---|
669 | false // preserve the opposite triangular part
|
---|
670 | >::run(m_matrix.const_cast_derived(), other.derived().nestedExpression());
|
---|
671 | }
|
---|
672 |
|
---|
673 | /***************************************************************************
|
---|
674 | * Implementation of TriangularBase methods
|
---|
675 | ***************************************************************************/
|
---|
676 |
|
---|
677 | /** Assigns a triangular or selfadjoint matrix to a dense matrix.
|
---|
678 | * If the matrix is triangular, the opposite part is set to zero. */
|
---|
679 | template<typename Derived>
|
---|
680 | template<typename DenseDerived>
|
---|
681 | void TriangularBase<Derived>::evalTo(MatrixBase<DenseDerived> &other) const
|
---|
682 | {
|
---|
683 | if(internal::traits<Derived>::Flags & EvalBeforeAssigningBit)
|
---|
684 | {
|
---|
685 | typename internal::plain_matrix_type<Derived>::type other_evaluated(rows(), cols());
|
---|
686 | evalToLazy(other_evaluated);
|
---|
687 | other.derived().swap(other_evaluated);
|
---|
688 | }
|
---|
689 | else
|
---|
690 | evalToLazy(other.derived());
|
---|
691 | }
|
---|
692 |
|
---|
693 | /** Assigns a triangular or selfadjoint matrix to a dense matrix.
|
---|
694 | * If the matrix is triangular, the opposite part is set to zero. */
|
---|
695 | template<typename Derived>
|
---|
696 | template<typename DenseDerived>
|
---|
697 | void TriangularBase<Derived>::evalToLazy(MatrixBase<DenseDerived> &other) const
|
---|
698 | {
|
---|
699 | enum {
|
---|
700 | unroll = DenseDerived::SizeAtCompileTime != Dynamic
|
---|
701 | && internal::traits<Derived>::CoeffReadCost != Dynamic
|
---|
702 | && DenseDerived::SizeAtCompileTime * internal::traits<Derived>::CoeffReadCost / 2
|
---|
703 | <= EIGEN_UNROLLING_LIMIT
|
---|
704 | };
|
---|
705 | other.derived().resize(this->rows(), this->cols());
|
---|
706 |
|
---|
707 | internal::triangular_assignment_selector
|
---|
708 | <DenseDerived, typename internal::traits<Derived>::MatrixTypeNestedCleaned, Derived::Mode,
|
---|
709 | unroll ? int(DenseDerived::SizeAtCompileTime) : Dynamic,
|
---|
710 | true // clear the opposite triangular part
|
---|
711 | >::run(other.derived(), derived().nestedExpression());
|
---|
712 | }
|
---|
713 |
|
---|
714 | /***************************************************************************
|
---|
715 | * Implementation of TriangularView methods
|
---|
716 | ***************************************************************************/
|
---|
717 |
|
---|
718 | /***************************************************************************
|
---|
719 | * Implementation of MatrixBase methods
|
---|
720 | ***************************************************************************/
|
---|
721 |
|
---|
722 | #ifdef EIGEN2_SUPPORT
|
---|
723 |
|
---|
724 | // implementation of part<>(), including the SelfAdjoint case.
|
---|
725 |
|
---|
726 | namespace internal {
|
---|
727 | template<typename MatrixType, unsigned int Mode>
|
---|
728 | struct eigen2_part_return_type
|
---|
729 | {
|
---|
730 | typedef TriangularView<MatrixType, Mode> type;
|
---|
731 | };
|
---|
732 |
|
---|
733 | template<typename MatrixType>
|
---|
734 | struct eigen2_part_return_type<MatrixType, SelfAdjoint>
|
---|
735 | {
|
---|
736 | typedef SelfAdjointView<MatrixType, Upper> type;
|
---|
737 | };
|
---|
738 | }
|
---|
739 |
|
---|
740 | /** \deprecated use MatrixBase::triangularView() */
|
---|
741 | template<typename Derived>
|
---|
742 | template<unsigned int Mode>
|
---|
743 | const typename internal::eigen2_part_return_type<Derived, Mode>::type MatrixBase<Derived>::part() const
|
---|
744 | {
|
---|
745 | return derived();
|
---|
746 | }
|
---|
747 |
|
---|
748 | /** \deprecated use MatrixBase::triangularView() */
|
---|
749 | template<typename Derived>
|
---|
750 | template<unsigned int Mode>
|
---|
751 | typename internal::eigen2_part_return_type<Derived, Mode>::type MatrixBase<Derived>::part()
|
---|
752 | {
|
---|
753 | return derived();
|
---|
754 | }
|
---|
755 | #endif
|
---|
756 |
|
---|
757 | /**
|
---|
758 | * \returns an expression of a triangular view extracted from the current matrix
|
---|
759 | *
|
---|
760 | * The parameter \a Mode can have the following values: \c #Upper, \c #StrictlyUpper, \c #UnitUpper,
|
---|
761 | * \c #Lower, \c #StrictlyLower, \c #UnitLower.
|
---|
762 | *
|
---|
763 | * Example: \include MatrixBase_extract.cpp
|
---|
764 | * Output: \verbinclude MatrixBase_extract.out
|
---|
765 | *
|
---|
766 | * \sa class TriangularView
|
---|
767 | */
|
---|
768 | template<typename Derived>
|
---|
769 | template<unsigned int Mode>
|
---|
770 | typename MatrixBase<Derived>::template TriangularViewReturnType<Mode>::Type
|
---|
771 | MatrixBase<Derived>::triangularView()
|
---|
772 | {
|
---|
773 | return derived();
|
---|
774 | }
|
---|
775 |
|
---|
776 | /** This is the const version of MatrixBase::triangularView() */
|
---|
777 | template<typename Derived>
|
---|
778 | template<unsigned int Mode>
|
---|
779 | typename MatrixBase<Derived>::template ConstTriangularViewReturnType<Mode>::Type
|
---|
780 | MatrixBase<Derived>::triangularView() const
|
---|
781 | {
|
---|
782 | return derived();
|
---|
783 | }
|
---|
784 |
|
---|
785 | /** \returns true if *this is approximately equal to an upper triangular matrix,
|
---|
786 | * within the precision given by \a prec.
|
---|
787 | *
|
---|
788 | * \sa isLowerTriangular()
|
---|
789 | */
|
---|
790 | template<typename Derived>
|
---|
791 | bool MatrixBase<Derived>::isUpperTriangular(const RealScalar& prec) const
|
---|
792 | {
|
---|
793 | using std::abs;
|
---|
794 | RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1);
|
---|
795 | for(Index j = 0; j < cols(); ++j)
|
---|
796 | {
|
---|
797 | Index maxi = (std::min)(j, rows()-1);
|
---|
798 | for(Index i = 0; i <= maxi; ++i)
|
---|
799 | {
|
---|
800 | RealScalar absValue = abs(coeff(i,j));
|
---|
801 | if(absValue > maxAbsOnUpperPart) maxAbsOnUpperPart = absValue;
|
---|
802 | }
|
---|
803 | }
|
---|
804 | RealScalar threshold = maxAbsOnUpperPart * prec;
|
---|
805 | for(Index j = 0; j < cols(); ++j)
|
---|
806 | for(Index i = j+1; i < rows(); ++i)
|
---|
807 | if(abs(coeff(i, j)) > threshold) return false;
|
---|
808 | return true;
|
---|
809 | }
|
---|
810 |
|
---|
811 | /** \returns true if *this is approximately equal to a lower triangular matrix,
|
---|
812 | * within the precision given by \a prec.
|
---|
813 | *
|
---|
814 | * \sa isUpperTriangular()
|
---|
815 | */
|
---|
816 | template<typename Derived>
|
---|
817 | bool MatrixBase<Derived>::isLowerTriangular(const RealScalar& prec) const
|
---|
818 | {
|
---|
819 | using std::abs;
|
---|
820 | RealScalar maxAbsOnLowerPart = static_cast<RealScalar>(-1);
|
---|
821 | for(Index j = 0; j < cols(); ++j)
|
---|
822 | for(Index i = j; i < rows(); ++i)
|
---|
823 | {
|
---|
824 | RealScalar absValue = abs(coeff(i,j));
|
---|
825 | if(absValue > maxAbsOnLowerPart) maxAbsOnLowerPart = absValue;
|
---|
826 | }
|
---|
827 | RealScalar threshold = maxAbsOnLowerPart * prec;
|
---|
828 | for(Index j = 1; j < cols(); ++j)
|
---|
829 | {
|
---|
830 | Index maxi = (std::min)(j, rows()-1);
|
---|
831 | for(Index i = 0; i < maxi; ++i)
|
---|
832 | if(abs(coeff(i, j)) > threshold) return false;
|
---|
833 | }
|
---|
834 | return true;
|
---|
835 | }
|
---|
836 |
|
---|
837 | } // end namespace Eigen
|
---|
838 |
|
---|
839 | #endif // EIGEN_TRIANGULARMATRIX_H
|
---|