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