DiagonalMatrix.h (14670B)
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 DiagonalVectorType::RealScalar RealScalar; 24 typedef typename internal::traits<Derived>::StorageKind StorageKind; 25 typedef typename internal::traits<Derived>::StorageIndex StorageIndex; 26 27 enum { 28 RowsAtCompileTime = DiagonalVectorType::SizeAtCompileTime, 29 ColsAtCompileTime = DiagonalVectorType::SizeAtCompileTime, 30 MaxRowsAtCompileTime = DiagonalVectorType::MaxSizeAtCompileTime, 31 MaxColsAtCompileTime = DiagonalVectorType::MaxSizeAtCompileTime, 32 IsVectorAtCompileTime = 0, 33 Flags = NoPreferredStorageOrderBit 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 EIGEN_DEVICE_FUNC 41 inline const Derived& derived() const { return *static_cast<const Derived*>(this); } 42 EIGEN_DEVICE_FUNC 43 inline Derived& derived() { return *static_cast<Derived*>(this); } 44 45 EIGEN_DEVICE_FUNC 46 DenseMatrixType toDenseMatrix() const { return derived(); } 47 48 EIGEN_DEVICE_FUNC 49 inline const DiagonalVectorType& diagonal() const { return derived().diagonal(); } 50 EIGEN_DEVICE_FUNC 51 inline DiagonalVectorType& diagonal() { return derived().diagonal(); } 52 53 EIGEN_DEVICE_FUNC 54 inline Index rows() const { return diagonal().size(); } 55 EIGEN_DEVICE_FUNC 56 inline Index cols() const { return diagonal().size(); } 57 58 template<typename MatrixDerived> 59 EIGEN_DEVICE_FUNC 60 const Product<Derived,MatrixDerived,LazyProduct> 61 operator*(const MatrixBase<MatrixDerived> &matrix) const 62 { 63 return Product<Derived, MatrixDerived, LazyProduct>(derived(),matrix.derived()); 64 } 65 66 typedef DiagonalWrapper<const CwiseUnaryOp<internal::scalar_inverse_op<Scalar>, const DiagonalVectorType> > InverseReturnType; 67 EIGEN_DEVICE_FUNC 68 inline const InverseReturnType 69 inverse() const 70 { 71 return InverseReturnType(diagonal().cwiseInverse()); 72 } 73 74 EIGEN_DEVICE_FUNC 75 inline const DiagonalWrapper<const EIGEN_EXPR_BINARYOP_SCALAR_RETURN_TYPE(DiagonalVectorType,Scalar,product) > 76 operator*(const Scalar& scalar) const 77 { 78 return DiagonalWrapper<const EIGEN_EXPR_BINARYOP_SCALAR_RETURN_TYPE(DiagonalVectorType,Scalar,product) >(diagonal() * scalar); 79 } 80 EIGEN_DEVICE_FUNC 81 friend inline const DiagonalWrapper<const EIGEN_SCALAR_BINARYOP_EXPR_RETURN_TYPE(Scalar,DiagonalVectorType,product) > 82 operator*(const Scalar& scalar, const DiagonalBase& other) 83 { 84 return DiagonalWrapper<const EIGEN_SCALAR_BINARYOP_EXPR_RETURN_TYPE(Scalar,DiagonalVectorType,product) >(scalar * other.diagonal()); 85 } 86 87 template<typename OtherDerived> 88 EIGEN_DEVICE_FUNC 89 #ifdef EIGEN_PARSED_BY_DOXYGEN 90 inline unspecified_expression_type 91 #else 92 inline const DiagonalWrapper<const EIGEN_CWISE_BINARY_RETURN_TYPE(DiagonalVectorType,typename OtherDerived::DiagonalVectorType,sum) > 93 #endif 94 operator+(const DiagonalBase<OtherDerived>& other) const 95 { 96 return (diagonal() + other.diagonal()).asDiagonal(); 97 } 98 99 template<typename OtherDerived> 100 EIGEN_DEVICE_FUNC 101 #ifdef EIGEN_PARSED_BY_DOXYGEN 102 inline unspecified_expression_type 103 #else 104 inline const DiagonalWrapper<const EIGEN_CWISE_BINARY_RETURN_TYPE(DiagonalVectorType,typename OtherDerived::DiagonalVectorType,difference) > 105 #endif 106 operator-(const DiagonalBase<OtherDerived>& other) const 107 { 108 return (diagonal() - other.diagonal()).asDiagonal(); 109 } 110 }; 111 112 #endif 113 114 /** \class DiagonalMatrix 115 * \ingroup Core_Module 116 * 117 * \brief Represents a diagonal matrix with its storage 118 * 119 * \param _Scalar the type of coefficients 120 * \param SizeAtCompileTime the dimension of the matrix, or Dynamic 121 * \param MaxSizeAtCompileTime the dimension of the matrix, or Dynamic. This parameter is optional and defaults 122 * to SizeAtCompileTime. Most of the time, you do not need to specify it. 123 * 124 * \sa class DiagonalWrapper 125 */ 126 127 namespace internal { 128 template<typename _Scalar, int SizeAtCompileTime, int MaxSizeAtCompileTime> 129 struct traits<DiagonalMatrix<_Scalar,SizeAtCompileTime,MaxSizeAtCompileTime> > 130 : traits<Matrix<_Scalar,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> > 131 { 132 typedef Matrix<_Scalar,SizeAtCompileTime,1,0,MaxSizeAtCompileTime,1> DiagonalVectorType; 133 typedef DiagonalShape StorageKind; 134 enum { 135 Flags = LvalueBit | NoPreferredStorageOrderBit 136 }; 137 }; 138 } 139 template<typename _Scalar, int SizeAtCompileTime, int MaxSizeAtCompileTime> 140 class DiagonalMatrix 141 : public DiagonalBase<DiagonalMatrix<_Scalar,SizeAtCompileTime,MaxSizeAtCompileTime> > 142 { 143 public: 144 #ifndef EIGEN_PARSED_BY_DOXYGEN 145 typedef typename internal::traits<DiagonalMatrix>::DiagonalVectorType DiagonalVectorType; 146 typedef const DiagonalMatrix& Nested; 147 typedef _Scalar Scalar; 148 typedef typename internal::traits<DiagonalMatrix>::StorageKind StorageKind; 149 typedef typename internal::traits<DiagonalMatrix>::StorageIndex StorageIndex; 150 #endif 151 152 protected: 153 154 DiagonalVectorType m_diagonal; 155 156 public: 157 158 /** const version of diagonal(). */ 159 EIGEN_DEVICE_FUNC 160 inline const DiagonalVectorType& diagonal() const { return m_diagonal; } 161 /** \returns a reference to the stored vector of diagonal coefficients. */ 162 EIGEN_DEVICE_FUNC 163 inline DiagonalVectorType& diagonal() { return m_diagonal; } 164 165 /** Default constructor without initialization */ 166 EIGEN_DEVICE_FUNC 167 inline DiagonalMatrix() {} 168 169 /** Constructs a diagonal matrix with given dimension */ 170 EIGEN_DEVICE_FUNC 171 explicit inline DiagonalMatrix(Index dim) : m_diagonal(dim) {} 172 173 /** 2D constructor. */ 174 EIGEN_DEVICE_FUNC 175 inline DiagonalMatrix(const Scalar& x, const Scalar& y) : m_diagonal(x,y) {} 176 177 /** 3D constructor. */ 178 EIGEN_DEVICE_FUNC 179 inline DiagonalMatrix(const Scalar& x, const Scalar& y, const Scalar& z) : m_diagonal(x,y,z) {} 180 181 #if EIGEN_HAS_CXX11 182 /** \brief Construct a diagonal matrix with fixed size from an arbitrary number of coefficients. \cpp11 183 * 184 * There exists C++98 anologue constructors for fixed-size diagonal matrices having 2 or 3 coefficients. 185 * 186 * \warning To construct a diagonal matrix of fixed size, the number of values passed to this 187 * constructor must match the fixed dimension of \c *this. 188 * 189 * \sa DiagonalMatrix(const Scalar&, const Scalar&) 190 * \sa DiagonalMatrix(const Scalar&, const Scalar&, const Scalar&) 191 */ 192 template <typename... ArgTypes> 193 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE 194 DiagonalMatrix(const Scalar& a0, const Scalar& a1, const Scalar& a2, const ArgTypes&... args) 195 : m_diagonal(a0, a1, a2, args...) {} 196 197 /** \brief Constructs a DiagonalMatrix and initializes it by elements given by an initializer list of initializer 198 * lists \cpp11 199 */ 200 EIGEN_DEVICE_FUNC 201 explicit EIGEN_STRONG_INLINE DiagonalMatrix(const std::initializer_list<std::initializer_list<Scalar>>& list) 202 : m_diagonal(list) {} 203 #endif // EIGEN_HAS_CXX11 204 205 /** Copy constructor. */ 206 template<typename OtherDerived> 207 EIGEN_DEVICE_FUNC 208 inline DiagonalMatrix(const DiagonalBase<OtherDerived>& other) : m_diagonal(other.diagonal()) {} 209 210 #ifndef EIGEN_PARSED_BY_DOXYGEN 211 /** copy constructor. prevent a default copy constructor from hiding the other templated constructor */ 212 inline DiagonalMatrix(const DiagonalMatrix& other) : m_diagonal(other.diagonal()) {} 213 #endif 214 215 /** generic constructor from expression of the diagonal coefficients */ 216 template<typename OtherDerived> 217 EIGEN_DEVICE_FUNC 218 explicit inline DiagonalMatrix(const MatrixBase<OtherDerived>& other) : m_diagonal(other) 219 {} 220 221 /** Copy operator. */ 222 template<typename OtherDerived> 223 EIGEN_DEVICE_FUNC 224 DiagonalMatrix& operator=(const DiagonalBase<OtherDerived>& other) 225 { 226 m_diagonal = other.diagonal(); 227 return *this; 228 } 229 230 #ifndef EIGEN_PARSED_BY_DOXYGEN 231 /** This is a special case of the templated operator=. Its purpose is to 232 * prevent a default operator= from hiding the templated operator=. 233 */ 234 EIGEN_DEVICE_FUNC 235 DiagonalMatrix& operator=(const DiagonalMatrix& other) 236 { 237 m_diagonal = other.diagonal(); 238 return *this; 239 } 240 #endif 241 242 /** Resizes to given size. */ 243 EIGEN_DEVICE_FUNC 244 inline void resize(Index size) { m_diagonal.resize(size); } 245 /** Sets all coefficients to zero. */ 246 EIGEN_DEVICE_FUNC 247 inline void setZero() { m_diagonal.setZero(); } 248 /** Resizes and sets all coefficients to zero. */ 249 EIGEN_DEVICE_FUNC 250 inline void setZero(Index size) { m_diagonal.setZero(size); } 251 /** Sets this matrix to be the identity matrix of the current size. */ 252 EIGEN_DEVICE_FUNC 253 inline void setIdentity() { m_diagonal.setOnes(); } 254 /** Sets this matrix to be the identity matrix of the given size. */ 255 EIGEN_DEVICE_FUNC 256 inline void setIdentity(Index size) { m_diagonal.setOnes(size); } 257 }; 258 259 /** \class DiagonalWrapper 260 * \ingroup Core_Module 261 * 262 * \brief Expression of a diagonal matrix 263 * 264 * \param _DiagonalVectorType the type of the vector of diagonal coefficients 265 * 266 * This class is an expression of a diagonal matrix, but not storing its own vector of diagonal coefficients, 267 * instead wrapping an existing vector expression. It is the return type of MatrixBase::asDiagonal() 268 * and most of the time this is the only way that it is used. 269 * 270 * \sa class DiagonalMatrix, class DiagonalBase, MatrixBase::asDiagonal() 271 */ 272 273 namespace internal { 274 template<typename _DiagonalVectorType> 275 struct traits<DiagonalWrapper<_DiagonalVectorType> > 276 { 277 typedef _DiagonalVectorType DiagonalVectorType; 278 typedef typename DiagonalVectorType::Scalar Scalar; 279 typedef typename DiagonalVectorType::StorageIndex StorageIndex; 280 typedef DiagonalShape StorageKind; 281 typedef typename traits<DiagonalVectorType>::XprKind XprKind; 282 enum { 283 RowsAtCompileTime = DiagonalVectorType::SizeAtCompileTime, 284 ColsAtCompileTime = DiagonalVectorType::SizeAtCompileTime, 285 MaxRowsAtCompileTime = DiagonalVectorType::MaxSizeAtCompileTime, 286 MaxColsAtCompileTime = DiagonalVectorType::MaxSizeAtCompileTime, 287 Flags = (traits<DiagonalVectorType>::Flags & LvalueBit) | NoPreferredStorageOrderBit 288 }; 289 }; 290 } 291 292 template<typename _DiagonalVectorType> 293 class DiagonalWrapper 294 : public DiagonalBase<DiagonalWrapper<_DiagonalVectorType> >, internal::no_assignment_operator 295 { 296 public: 297 #ifndef EIGEN_PARSED_BY_DOXYGEN 298 typedef _DiagonalVectorType DiagonalVectorType; 299 typedef DiagonalWrapper Nested; 300 #endif 301 302 /** Constructor from expression of diagonal coefficients to wrap. */ 303 EIGEN_DEVICE_FUNC 304 explicit inline DiagonalWrapper(DiagonalVectorType& a_diagonal) : m_diagonal(a_diagonal) {} 305 306 /** \returns a const reference to the wrapped expression of diagonal coefficients. */ 307 EIGEN_DEVICE_FUNC 308 const DiagonalVectorType& diagonal() const { return m_diagonal; } 309 310 protected: 311 typename DiagonalVectorType::Nested m_diagonal; 312 }; 313 314 /** \returns a pseudo-expression of a diagonal matrix with *this as vector of diagonal coefficients 315 * 316 * \only_for_vectors 317 * 318 * Example: \include MatrixBase_asDiagonal.cpp 319 * Output: \verbinclude MatrixBase_asDiagonal.out 320 * 321 * \sa class DiagonalWrapper, class DiagonalMatrix, diagonal(), isDiagonal() 322 **/ 323 template<typename Derived> 324 EIGEN_DEVICE_FUNC inline const DiagonalWrapper<const Derived> 325 MatrixBase<Derived>::asDiagonal() const 326 { 327 return DiagonalWrapper<const Derived>(derived()); 328 } 329 330 /** \returns true if *this is approximately equal to a diagonal matrix, 331 * within the precision given by \a prec. 332 * 333 * Example: \include MatrixBase_isDiagonal.cpp 334 * Output: \verbinclude MatrixBase_isDiagonal.out 335 * 336 * \sa asDiagonal() 337 */ 338 template<typename Derived> 339 bool MatrixBase<Derived>::isDiagonal(const RealScalar& prec) const 340 { 341 if(cols() != rows()) return false; 342 RealScalar maxAbsOnDiagonal = static_cast<RealScalar>(-1); 343 for(Index j = 0; j < cols(); ++j) 344 { 345 RealScalar absOnDiagonal = numext::abs(coeff(j,j)); 346 if(absOnDiagonal > maxAbsOnDiagonal) maxAbsOnDiagonal = absOnDiagonal; 347 } 348 for(Index j = 0; j < cols(); ++j) 349 for(Index i = 0; i < j; ++i) 350 { 351 if(!internal::isMuchSmallerThan(coeff(i, j), maxAbsOnDiagonal, prec)) return false; 352 if(!internal::isMuchSmallerThan(coeff(j, i), maxAbsOnDiagonal, prec)) return false; 353 } 354 return true; 355 } 356 357 namespace internal { 358 359 template<> struct storage_kind_to_shape<DiagonalShape> { typedef DiagonalShape Shape; }; 360 361 struct Diagonal2Dense {}; 362 363 template<> struct AssignmentKind<DenseShape,DiagonalShape> { typedef Diagonal2Dense Kind; }; 364 365 // Diagonal matrix to Dense assignment 366 template< typename DstXprType, typename SrcXprType, typename Functor> 367 struct Assignment<DstXprType, SrcXprType, Functor, Diagonal2Dense> 368 { 369 static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &/*func*/) 370 { 371 Index dstRows = src.rows(); 372 Index dstCols = src.cols(); 373 if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) 374 dst.resize(dstRows, dstCols); 375 376 dst.setZero(); 377 dst.diagonal() = src.diagonal(); 378 } 379 380 static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &/*func*/) 381 { dst.diagonal() += src.diagonal(); } 382 383 static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &/*func*/) 384 { dst.diagonal() -= src.diagonal(); } 385 }; 386 387 } // namespace internal 388 389 } // end namespace Eigen 390 391 #endif // EIGEN_DIAGONALMATRIX_H