Yoji KURODA / Eigen

Eigne Matrix Class Library

Dependents:   Eigen_test Odometry_test AttitudeEstimation_usingTicker MPU9250_Quaternion_Binary_Serial ... more

Eigen Matrix Class Library for mbed.

Finally, you can use Eigen on your mbed!!!

src/Core/TriangularMatrix.h@0:13a5d365ba16, 2016-10-13 (annotated)

Committer:
ykuroda
Date:
Thu Oct 13 04:07:23 2016 +0000
Revision:
0:13a5d365ba16
First commint, Eigne Matrix Class Library

Who changed what in which revision?

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