1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> 5 // Copyright (C) 2007-2009 Benoit Jacob <jacob.benoit.1@gmail.com> 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_DIAGONALMATRIX_H 12 #define EIGEN_DIAGONALMATRIX_H 13 14 namespace Eigen { 15 16 #ifndef EIGEN_PARSED_BY_DOXYGEN 17 template<typename Derived> 18 class DiagonalBase : public EigenBase<Derived> 19 { 20 public: 21 typedef typename internal::traits<Derived>::DiagonalVectorType DiagonalVectorType; 22 typedef typename DiagonalVectorType::Scalar Scalar; 23 typedef typename internal::traits<Derived>::StorageKind StorageKind; 24 typedef typename internal::traits<Derived>::Index Index; 25 26 enum { 27 RowsAtCompileTime = DiagonalVectorType::SizeAtCompileTime, 28 ColsAtCompileTime = DiagonalVectorType::SizeAtCompileTime, 29 MaxRowsAtCompileTime = DiagonalVectorType::MaxSizeAtCompileTime, 30 MaxColsAtCompileTime = DiagonalVectorType::MaxSizeAtCompileTime, 31 IsVectorAtCompileTime = 0, 32 Flags = 0 33 }; 34 35 typedef Matrix<Scalar, RowsAtCompileTime, ColsAtCompileTime, 0, MaxRowsAtCompileTime, MaxColsAtCompileTime> DenseMatrixType; 36 typedef DenseMatrixType DenseType; 37 typedef DiagonalMatrix<Scalar,DiagonalVectorType::SizeAtCompileTime,DiagonalVectorType::MaxSizeAtCompileTime> PlainObject; 38 derived()39 inline const Derived& derived() const { return *static_cast<const Derived*>(this); } derived()40 inline Derived& derived() { return *static_cast<Derived*>(this); } 41 toDenseMatrix()42 DenseMatrixType toDenseMatrix() const { return derived(); } 43 template<typename DenseDerived> 44 void evalTo(MatrixBase<DenseDerived> &other) const; 45 template<typename DenseDerived> addTo(MatrixBase<DenseDerived> & other)46 void addTo(MatrixBase<DenseDerived> &other) const 47 { other.diagonal() += diagonal(); } 48 template<typename DenseDerived> subTo(MatrixBase<DenseDerived> & other)49 void subTo(MatrixBase<DenseDerived> &other) const 50 { other.diagonal() -= diagonal(); } 51 diagonal()52 inline const DiagonalVectorType& diagonal() const { return derived().diagonal(); } diagonal()53 inline DiagonalVectorType& diagonal() { return derived().diagonal(); } 54 rows()55 inline Index rows() const { return diagonal().size(); } cols()56 inline Index cols() const { return diagonal().size(); } 57 58 template<typename MatrixDerived> 59 const DiagonalProduct<MatrixDerived, Derived, OnTheLeft> 60 operator*(const MatrixBase<MatrixDerived> &matrix) const; 61 62 inline const DiagonalWrapper<const CwiseUnaryOp<internal::scalar_inverse_op<Scalar>, const DiagonalVectorType> > inverse()63 inverse() const 64 { 65 return diagonal().cwiseInverse(); 66 } 67 68 #ifdef EIGEN2_SUPPORT 69 template<typename OtherDerived> 70 bool isApprox(const DiagonalBase<OtherDerived>& other, typename NumTraits<Scalar>::Real precision = NumTraits<Scalar>::dummy_precision()) const 71 { 72 return diagonal().isApprox(other.diagonal(), precision); 73 } 74 template<typename OtherDerived> 75 bool isApprox(const MatrixBase<OtherDerived>& other, typename NumTraits<Scalar>::Real precision = NumTraits<Scalar>::dummy_precision()) const 76 { 77 return toDenseMatrix().isApprox(other, precision); 78 } 79 #endif 80 }; 81 82 template<typename Derived> 83 template<typename DenseDerived> evalTo(MatrixBase<DenseDerived> & other)84 void DiagonalBase<Derived>::evalTo(MatrixBase<DenseDerived> &other) const 85 { 86 other.setZero(); 87 other.diagonal() = diagonal(); 88 } 89 #endif 90 91 /** \class DiagonalMatrix 92 * \ingroup Core_Module 93 * 94 * \brief Represents a diagonal matrix with its storage 95 * 96 * \param _Scalar the type of coefficients 97 * \param SizeAtCompileTime the dimension of the matrix, or Dynamic 98 * \param MaxSizeAtCompileTime the dimension of the matrix, or Dynamic. This parameter is optional and defaults 99 * to SizeAtCompileTime. Most of the time, you do not need to specify it. 100 * 101 * \sa class DiagonalWrapper 102 */ 103 104 namespace internal { 105 template<typename _Scalar, int SizeAtCompileTime, int MaxSizeAtCompileTime> 106 struct traits<DiagonalMatrix<_Scalar,SizeAtCompileTime,MaxSizeAtCompileTime> > 107 : traits<Matrix<_Scalar,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> > 108 { 109 typedef Matrix<_Scalar,SizeAtCompileTime,1,0,MaxSizeAtCompileTime,1> DiagonalVectorType; 110 typedef Dense StorageKind; 111 typedef DenseIndex Index; 112 enum { 113 Flags = LvalueBit 114 }; 115 }; 116 } 117 template<typename _Scalar, int SizeAtCompileTime, int MaxSizeAtCompileTime> 118 class DiagonalMatrix 119 : public DiagonalBase<DiagonalMatrix<_Scalar,SizeAtCompileTime,MaxSizeAtCompileTime> > 120 { 121 public: 122 #ifndef EIGEN_PARSED_BY_DOXYGEN 123 typedef typename internal::traits<DiagonalMatrix>::DiagonalVectorType DiagonalVectorType; 124 typedef const DiagonalMatrix& Nested; 125 typedef _Scalar Scalar; 126 typedef typename internal::traits<DiagonalMatrix>::StorageKind StorageKind; 127 typedef typename internal::traits<DiagonalMatrix>::Index Index; 128 #endif 129 130 protected: 131 132 DiagonalVectorType m_diagonal; 133 134 public: 135 136 /** const version of diagonal(). */ 137 inline const DiagonalVectorType& diagonal() const { return m_diagonal; } 138 /** \returns a reference to the stored vector of diagonal coefficients. */ 139 inline DiagonalVectorType& diagonal() { return m_diagonal; } 140 141 /** Default constructor without initialization */ 142 inline DiagonalMatrix() {} 143 144 /** Constructs a diagonal matrix with given dimension */ 145 inline DiagonalMatrix(Index dim) : m_diagonal(dim) {} 146 147 /** 2D constructor. */ 148 inline DiagonalMatrix(const Scalar& x, const Scalar& y) : m_diagonal(x,y) {} 149 150 /** 3D constructor. */ 151 inline DiagonalMatrix(const Scalar& x, const Scalar& y, const Scalar& z) : m_diagonal(x,y,z) {} 152 153 /** Copy constructor. */ 154 template<typename OtherDerived> 155 inline DiagonalMatrix(const DiagonalBase<OtherDerived>& other) : m_diagonal(other.diagonal()) {} 156 157 #ifndef EIGEN_PARSED_BY_DOXYGEN 158 /** copy constructor. prevent a default copy constructor from hiding the other templated constructor */ 159 inline DiagonalMatrix(const DiagonalMatrix& other) : m_diagonal(other.diagonal()) {} 160 #endif 161 162 /** generic constructor from expression of the diagonal coefficients */ 163 template<typename OtherDerived> 164 explicit inline DiagonalMatrix(const MatrixBase<OtherDerived>& other) : m_diagonal(other) 165 {} 166 167 /** Copy operator. */ 168 template<typename OtherDerived> 169 DiagonalMatrix& operator=(const DiagonalBase<OtherDerived>& other) 170 { 171 m_diagonal = other.diagonal(); 172 return *this; 173 } 174 175 #ifndef EIGEN_PARSED_BY_DOXYGEN 176 /** This is a special case of the templated operator=. Its purpose is to 177 * prevent a default operator= from hiding the templated operator=. 178 */ 179 DiagonalMatrix& operator=(const DiagonalMatrix& other) 180 { 181 m_diagonal = other.diagonal(); 182 return *this; 183 } 184 #endif 185 186 /** Resizes to given size. */ 187 inline void resize(Index size) { m_diagonal.resize(size); } 188 /** Sets all coefficients to zero. */ 189 inline void setZero() { m_diagonal.setZero(); } 190 /** Resizes and sets all coefficients to zero. */ 191 inline void setZero(Index size) { m_diagonal.setZero(size); } 192 /** Sets this matrix to be the identity matrix of the current size. */ 193 inline void setIdentity() { m_diagonal.setOnes(); } 194 /** Sets this matrix to be the identity matrix of the given size. */ 195 inline void setIdentity(Index size) { m_diagonal.setOnes(size); } 196 }; 197 198 /** \class DiagonalWrapper 199 * \ingroup Core_Module 200 * 201 * \brief Expression of a diagonal matrix 202 * 203 * \param _DiagonalVectorType the type of the vector of diagonal coefficients 204 * 205 * This class is an expression of a diagonal matrix, but not storing its own vector of diagonal coefficients, 206 * instead wrapping an existing vector expression. It is the return type of MatrixBase::asDiagonal() 207 * and most of the time this is the only way that it is used. 208 * 209 * \sa class DiagonalMatrix, class DiagonalBase, MatrixBase::asDiagonal() 210 */ 211 212 namespace internal { 213 template<typename _DiagonalVectorType> 214 struct traits<DiagonalWrapper<_DiagonalVectorType> > 215 { 216 typedef _DiagonalVectorType DiagonalVectorType; 217 typedef typename DiagonalVectorType::Scalar Scalar; 218 typedef typename DiagonalVectorType::Index Index; 219 typedef typename DiagonalVectorType::StorageKind StorageKind; 220 enum { 221 RowsAtCompileTime = DiagonalVectorType::SizeAtCompileTime, 222 ColsAtCompileTime = DiagonalVectorType::SizeAtCompileTime, 223 MaxRowsAtCompileTime = DiagonalVectorType::SizeAtCompileTime, 224 MaxColsAtCompileTime = DiagonalVectorType::SizeAtCompileTime, 225 Flags = traits<DiagonalVectorType>::Flags & LvalueBit 226 }; 227 }; 228 } 229 230 template<typename _DiagonalVectorType> 231 class DiagonalWrapper 232 : public DiagonalBase<DiagonalWrapper<_DiagonalVectorType> >, internal::no_assignment_operator 233 { 234 public: 235 #ifndef EIGEN_PARSED_BY_DOXYGEN 236 typedef _DiagonalVectorType DiagonalVectorType; 237 typedef DiagonalWrapper Nested; 238 #endif 239 240 /** Constructor from expression of diagonal coefficients to wrap. */ 241 inline DiagonalWrapper(DiagonalVectorType& diagonal) : m_diagonal(diagonal) {} 242 243 /** \returns a const reference to the wrapped expression of diagonal coefficients. */ 244 const DiagonalVectorType& diagonal() const { return m_diagonal; } 245 246 protected: 247 typename DiagonalVectorType::Nested m_diagonal; 248 }; 249 250 /** \returns a pseudo-expression of a diagonal matrix with *this as vector of diagonal coefficients 251 * 252 * \only_for_vectors 253 * 254 * Example: \include MatrixBase_asDiagonal.cpp 255 * Output: \verbinclude MatrixBase_asDiagonal.out 256 * 257 * \sa class DiagonalWrapper, class DiagonalMatrix, diagonal(), isDiagonal() 258 **/ 259 template<typename Derived> 260 inline const DiagonalWrapper<const Derived> 261 MatrixBase<Derived>::asDiagonal() const 262 { 263 return derived(); 264 } 265 266 /** \returns true if *this is approximately equal to a diagonal matrix, 267 * within the precision given by \a prec. 268 * 269 * Example: \include MatrixBase_isDiagonal.cpp 270 * Output: \verbinclude MatrixBase_isDiagonal.out 271 * 272 * \sa asDiagonal() 273 */ 274 template<typename Derived> 275 bool MatrixBase<Derived>::isDiagonal(RealScalar prec) const 276 { 277 if(cols() != rows()) return false; 278 RealScalar maxAbsOnDiagonal = static_cast<RealScalar>(-1); 279 for(Index j = 0; j < cols(); ++j) 280 { 281 RealScalar absOnDiagonal = internal::abs(coeff(j,j)); 282 if(absOnDiagonal > maxAbsOnDiagonal) maxAbsOnDiagonal = absOnDiagonal; 283 } 284 for(Index j = 0; j < cols(); ++j) 285 for(Index i = 0; i < j; ++i) 286 { 287 if(!internal::isMuchSmallerThan(coeff(i, j), maxAbsOnDiagonal, prec)) return false; 288 if(!internal::isMuchSmallerThan(coeff(j, i), maxAbsOnDiagonal, prec)) return false; 289 } 290 return true; 291 } 292 293 } // end namespace Eigen 294 295 #endif // EIGEN_DIAGONALMATRIX_H 296