src/Core/DiagonalMatrix.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 // Copyright (C) 2007-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
0006 //
0007 // This Source Code Form is subject to the terms of the Mozilla
0008 // Public License v. 2.0. If a copy of the MPL was not distributed
0009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0010
0011 #ifndef EIGEN_DIAGONALMATRIX_H
0012 #define EIGEN_DIAGONALMATRIX_H
0013
0014 namespace Eigen {
0015
0016 #ifndef EIGEN_PARSED_BY_DOXYGEN
0017 template<typename Derived>
0018 class DiagonalBase : public EigenBase<Derived>
0019 {
0020   public:
0021     typedef typename internal::traits<Derived>::DiagonalVectorType DiagonalVectorType;
0022     typedef typename DiagonalVectorType::Scalar Scalar;
0023     typedef typename DiagonalVectorType::RealScalar RealScalar;
0024     typedef typename internal::traits<Derived>::StorageKind StorageKind;
0025     typedef typename internal::traits<Derived>::StorageIndex StorageIndex;
0026
0027     enum {
0028       RowsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
0029       ColsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
0030       MaxRowsAtCompileTime = DiagonalVectorType::MaxSizeAtCompileTime,
0031       MaxColsAtCompileTime = DiagonalVectorType::MaxSizeAtCompileTime,
0032       IsVectorAtCompileTime = 0,
0033       Flags = NoPreferredStorageOrderBit
0034     };
0035
0036     typedef Matrix<Scalar, RowsAtCompileTime, ColsAtCompileTime, 0, MaxRowsAtCompileTime, MaxColsAtCompileTime> DenseMatrixType;
0037     typedef DenseMatrixType DenseType;
0038     typedef DiagonalMatrix<Scalar,DiagonalVectorType::SizeAtCompileTime,DiagonalVectorType::MaxSizeAtCompileTime> PlainObject;
0039
0040     EIGEN_DEVICE_FUNC
0041     inline const Derived& derived() const { return *static_cast<const Derived*>(this); }
0042     EIGEN_DEVICE_FUNC
0043     inline Derived& derived() { return *static_cast<Derived*>(this); }
0044
0045     EIGEN_DEVICE_FUNC
0046     DenseMatrixType toDenseMatrix() const { return derived(); }
0047
0048     EIGEN_DEVICE_FUNC
0049     inline const DiagonalVectorType& diagonal() const { return derived().diagonal(); }
0050     EIGEN_DEVICE_FUNC
0051     inline DiagonalVectorType& diagonal() { return derived().diagonal(); }
0052
0053     EIGEN_DEVICE_FUNC
0054     inline Index rows() const { return diagonal().size(); }
0055     EIGEN_DEVICE_FUNC
0056     inline Index cols() const { return diagonal().size(); }
0057
0058     template<typename MatrixDerived>
0059     EIGEN_DEVICE_FUNC
0060     const Product<Derived,MatrixDerived,LazyProduct>
0061     operator*(const MatrixBase<MatrixDerived> &matrix) const
0062     {
0063       return Product<Derived, MatrixDerived, LazyProduct>(derived(),matrix.derived());
0064     }
0065
0066     typedef DiagonalWrapper<const CwiseUnaryOp<internal::scalar_inverse_op<Scalar>, const DiagonalVectorType> > InverseReturnType;
0067     EIGEN_DEVICE_FUNC
0068     inline const InverseReturnType
0069     inverse() const
0070     {
0071       return InverseReturnType(diagonal().cwiseInverse());
0072     }
0073
0074     EIGEN_DEVICE_FUNC
0075     inline const DiagonalWrapper<const EIGEN_EXPR_BINARYOP_SCALAR_RETURN_TYPE(DiagonalVectorType,Scalar,product) >
0076     operator*(const Scalar& scalar) const
0077     {
0078       return DiagonalWrapper<const EIGEN_EXPR_BINARYOP_SCALAR_RETURN_TYPE(DiagonalVectorType,Scalar,product) >(diagonal() * scalar);
0079     }
0080     EIGEN_DEVICE_FUNC
0081     friend inline const DiagonalWrapper<const EIGEN_SCALAR_BINARYOP_EXPR_RETURN_TYPE(Scalar,DiagonalVectorType,product) >
0082     operator*(const Scalar& scalar, const DiagonalBase& other)
0083     {
0084       return DiagonalWrapper<const EIGEN_SCALAR_BINARYOP_EXPR_RETURN_TYPE(Scalar,DiagonalVectorType,product) >(scalar * other.diagonal());
0085     }
0086
0087     template<typename OtherDerived>
0088     EIGEN_DEVICE_FUNC
0089     #ifdef EIGEN_PARSED_BY_DOXYGEN
0090     inline unspecified_expression_type
0091     #else
0092     inline const DiagonalWrapper<const EIGEN_CWISE_BINARY_RETURN_TYPE(DiagonalVectorType,typename OtherDerived::DiagonalVectorType,sum) >
0093     #endif
0094     operator+(const DiagonalBase<OtherDerived>& other) const
0095     {
0096       return (diagonal() + other.diagonal()).asDiagonal();
0097     }
0098
0099     template<typename OtherDerived>
0100     EIGEN_DEVICE_FUNC
0101     #ifdef EIGEN_PARSED_BY_DOXYGEN
0102     inline unspecified_expression_type
0103     #else
0104     inline const DiagonalWrapper<const EIGEN_CWISE_BINARY_RETURN_TYPE(DiagonalVectorType,typename OtherDerived::DiagonalVectorType,difference) >
0105     #endif
0106     operator-(const DiagonalBase<OtherDerived>& other) const
0107     {
0108       return (diagonal() - other.diagonal()).asDiagonal();
0109     }
0110 };
0111
0112 #endif
0113
0114 /** \class DiagonalMatrix
0115   * \ingroup Core_Module
0116   *
0117   * \brief Represents a diagonal matrix with its storage
0118   *
0119   * \param _Scalar the type of coefficients
0120   * \param SizeAtCompileTime the dimension of the matrix, or Dynamic
0121   * \param MaxSizeAtCompileTime the dimension of the matrix, or Dynamic. This parameter is optional and defaults
0122   *        to SizeAtCompileTime. Most of the time, you do not need to specify it.
0123   *
0124   * \sa class DiagonalWrapper
0125   */
0126
0127 namespace internal {
0128 template<typename _Scalar, int SizeAtCompileTime, int MaxSizeAtCompileTime>
0129 struct traits<DiagonalMatrix<_Scalar,SizeAtCompileTime,MaxSizeAtCompileTime> >
0130  : traits<Matrix<_Scalar,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
0131 {
0132   typedef Matrix<_Scalar,SizeAtCompileTime,1,0,MaxSizeAtCompileTime,1> DiagonalVectorType;
0133   typedef DiagonalShape StorageKind;
0134   enum {
0135     Flags = LvalueBit | NoPreferredStorageOrderBit
0136   };
0137 };
0138 }
0139 template<typename _Scalar, int SizeAtCompileTime, int MaxSizeAtCompileTime>
0140 class DiagonalMatrix
0141   : public DiagonalBase<DiagonalMatrix<_Scalar,SizeAtCompileTime,MaxSizeAtCompileTime> >
0142 {
0143   public:
0144     #ifndef EIGEN_PARSED_BY_DOXYGEN
0145     typedef typename internal::traits<DiagonalMatrix>::DiagonalVectorType DiagonalVectorType;
0146     typedef const DiagonalMatrix& Nested;
0147     typedef _Scalar Scalar;
0148     typedef typename internal::traits<DiagonalMatrix>::StorageKind StorageKind;
0149     typedef typename internal::traits<DiagonalMatrix>::StorageIndex StorageIndex;
0150     #endif
0151
0152   protected:
0153
0154     DiagonalVectorType m_diagonal;
0155
0156   public:
0157
0158     /** const version of diagonal(). */
0159     EIGEN_DEVICE_FUNC
0160     inline const DiagonalVectorType& diagonal() const { return m_diagonal; }
0161     /** \returns a reference to the stored vector of diagonal coefficients. */
0162     EIGEN_DEVICE_FUNC
0163     inline DiagonalVectorType& diagonal() { return m_diagonal; }
0164
0165     /** Default constructor without initialization */
0166     EIGEN_DEVICE_FUNC
0167     inline DiagonalMatrix() {}
0168
0169     /** Constructs a diagonal matrix with given dimension  */
0170     EIGEN_DEVICE_FUNC
0171     explicit inline DiagonalMatrix(Index dim) : m_diagonal(dim) {}
0172
0173     /** 2D constructor. */
0174     EIGEN_DEVICE_FUNC
0175     inline DiagonalMatrix(const Scalar& x, const Scalar& y) : m_diagonal(x,y) {}
0176
0177     /** 3D constructor. */
0178     EIGEN_DEVICE_FUNC
0179     inline DiagonalMatrix(const Scalar& x, const Scalar& y, const Scalar& z) : m_diagonal(x,y,z) {}
0180
0181     #if EIGEN_HAS_CXX11
0182     /** \brief Construct a diagonal matrix with fixed size from an arbitrary number of coefficients. \cpp11
0183       *
0184       * There exists C++98 anologue constructors for fixed-size diagonal matrices having 2 or 3 coefficients.
0185       *
0186       * \warning To construct a diagonal matrix of fixed size, the number of values passed to this
0187       * constructor must match the fixed dimension of \c *this.
0188       *
0189       * \sa DiagonalMatrix(const Scalar&, const Scalar&)
0190       * \sa DiagonalMatrix(const Scalar&, const Scalar&, const Scalar&)
0191       */
0192     template <typename... ArgTypes>
0193     EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0194     DiagonalMatrix(const Scalar& a0, const Scalar& a1, const Scalar& a2, const ArgTypes&... args)
0195       : m_diagonal(a0, a1, a2, args...) {}
0196
0197     /** \brief Constructs a DiagonalMatrix and initializes it by elements given by an initializer list of initializer
0198       * lists \cpp11
0199       */
0200     EIGEN_DEVICE_FUNC
0201     explicit EIGEN_STRONG_INLINE DiagonalMatrix(const std::initializer_list<std::initializer_list<Scalar>>& list)
0202       : m_diagonal(list) {}
0203     #endif  // EIGEN_HAS_CXX11
0204
0205     /** Copy constructor. */
0206     template<typename OtherDerived>
0207     EIGEN_DEVICE_FUNC
0208     inline DiagonalMatrix(const DiagonalBase<OtherDerived>& other) : m_diagonal(other.diagonal()) {}
0209
0210     #ifndef EIGEN_PARSED_BY_DOXYGEN
0211     /** copy constructor. prevent a default copy constructor from hiding the other templated constructor */
0212     inline DiagonalMatrix(const DiagonalMatrix& other) : m_diagonal(other.diagonal()) {}
0213     #endif
0214
0215     /** generic constructor from expression of the diagonal coefficients */
0216     template<typename OtherDerived>
0217     EIGEN_DEVICE_FUNC
0218     explicit inline DiagonalMatrix(const MatrixBase<OtherDerived>& other) : m_diagonal(other)
0219     {}
0220
0221     /** Copy operator. */
0222     template<typename OtherDerived>
0223     EIGEN_DEVICE_FUNC
0224     DiagonalMatrix& operator=(const DiagonalBase<OtherDerived>& other)
0225     {
0226       m_diagonal = other.diagonal();
0227       return *this;
0228     }
0229
0230     #ifndef EIGEN_PARSED_BY_DOXYGEN
0231     /** This is a special case of the templated operator=. Its purpose is to
0232       * prevent a default operator= from hiding the templated operator=.
0233       */
0234     EIGEN_DEVICE_FUNC
0235     DiagonalMatrix& operator=(const DiagonalMatrix& other)
0236     {
0237       m_diagonal = other.diagonal();
0238       return *this;
0239     }
0240     #endif
0241
0242     /** Resizes to given size. */
0243     EIGEN_DEVICE_FUNC
0244     inline void resize(Index size) { m_diagonal.resize(size); }
0245     /** Sets all coefficients to zero. */
0246     EIGEN_DEVICE_FUNC
0247     inline void setZero() { m_diagonal.setZero(); }
0248     /** Resizes and sets all coefficients to zero. */
0249     EIGEN_DEVICE_FUNC
0250     inline void setZero(Index size) { m_diagonal.setZero(size); }
0251     /** Sets this matrix to be the identity matrix of the current size. */
0252     EIGEN_DEVICE_FUNC
0253     inline void setIdentity() { m_diagonal.setOnes(); }
0254     /** Sets this matrix to be the identity matrix of the given size. */
0255     EIGEN_DEVICE_FUNC
0256     inline void setIdentity(Index size) { m_diagonal.setOnes(size); }
0257 };
0258
0259 /** \class DiagonalWrapper
0260   * \ingroup Core_Module
0261   *
0262   * \brief Expression of a diagonal matrix
0263   *
0264   * \param _DiagonalVectorType the type of the vector of diagonal coefficients
0265   *
0266   * This class is an expression of a diagonal matrix, but not storing its own vector of diagonal coefficients,
0267   * instead wrapping an existing vector expression. It is the return type of MatrixBase::asDiagonal()
0268   * and most of the time this is the only way that it is used.
0269   *
0270   * \sa class DiagonalMatrix, class DiagonalBase, MatrixBase::asDiagonal()
0271   */
0272
0273 namespace internal {
0274 template<typename _DiagonalVectorType>
0275 struct traits<DiagonalWrapper<_DiagonalVectorType> >
0276 {
0277   typedef _DiagonalVectorType DiagonalVectorType;
0278   typedef typename DiagonalVectorType::Scalar Scalar;
0279   typedef typename DiagonalVectorType::StorageIndex StorageIndex;
0280   typedef DiagonalShape StorageKind;
0281   typedef typename traits<DiagonalVectorType>::XprKind XprKind;
0282   enum {
0283     RowsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
0284     ColsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
0285     MaxRowsAtCompileTime = DiagonalVectorType::MaxSizeAtCompileTime,
0286     MaxColsAtCompileTime = DiagonalVectorType::MaxSizeAtCompileTime,
0287     Flags =  (traits<DiagonalVectorType>::Flags & LvalueBit) | NoPreferredStorageOrderBit
0288   };
0289 };
0290 }
0291
0292 template<typename _DiagonalVectorType>
0293 class DiagonalWrapper
0294   : public DiagonalBase<DiagonalWrapper<_DiagonalVectorType> >, internal::no_assignment_operator
0295 {
0296   public:
0297     #ifndef EIGEN_PARSED_BY_DOXYGEN
0298     typedef _DiagonalVectorType DiagonalVectorType;
0299     typedef DiagonalWrapper Nested;
0300     #endif
0301
0302     /** Constructor from expression of diagonal coefficients to wrap. */
0303     EIGEN_DEVICE_FUNC
0304     explicit inline DiagonalWrapper(DiagonalVectorType& a_diagonal) : m_diagonal(a_diagonal) {}
0305
0306     /** \returns a const reference to the wrapped expression of diagonal coefficients. */
0307     EIGEN_DEVICE_FUNC
0308     const DiagonalVectorType& diagonal() const { return m_diagonal; }
0309
0310   protected:
0311     typename DiagonalVectorType::Nested m_diagonal;
0312 };
0313
0314 /** \returns a pseudo-expression of a diagonal matrix with *this as vector of diagonal coefficients
0315   *
0316   * \only_for_vectors
0317   *
0318   * Example: \include MatrixBase_asDiagonal.cpp
0319   * Output: \verbinclude MatrixBase_asDiagonal.out
0320   *
0321   * \sa class DiagonalWrapper, class DiagonalMatrix, diagonal(), isDiagonal()
0322   **/
0323 template<typename Derived>
0324 EIGEN_DEVICE_FUNC inline const DiagonalWrapper<const Derived>
0325 MatrixBase<Derived>::asDiagonal() const
0326 {
0327   return DiagonalWrapper<const Derived>(derived());
0328 }
0329
0330 /** \returns true if *this is approximately equal to a diagonal matrix,
0331   *          within the precision given by \a prec.
0332   *
0333   * Example: \include MatrixBase_isDiagonal.cpp
0334   * Output: \verbinclude MatrixBase_isDiagonal.out
0335   *
0336   * \sa asDiagonal()
0337   */
0338 template<typename Derived>
0339 bool MatrixBase<Derived>::isDiagonal(const RealScalar& prec) const
0340 {
0341   if(cols() != rows()) return false;
0342   RealScalar maxAbsOnDiagonal = static_cast<RealScalar>(-1);
0343   for(Index j = 0; j < cols(); ++j)
0344   {
0345     RealScalar absOnDiagonal = numext::abs(coeff(j,j));
0346     if(absOnDiagonal > maxAbsOnDiagonal) maxAbsOnDiagonal = absOnDiagonal;
0347   }
0348   for(Index j = 0; j < cols(); ++j)
0349     for(Index i = 0; i < j; ++i)
0350     {
0351       if(!internal::isMuchSmallerThan(coeff(i, j), maxAbsOnDiagonal, prec)) return false;
0352       if(!internal::isMuchSmallerThan(coeff(j, i), maxAbsOnDiagonal, prec)) return false;
0353     }
0354   return true;
0355 }
0356
0357 namespace internal {
0358
0359 template<> struct storage_kind_to_shape<DiagonalShape> { typedef DiagonalShape Shape; };
0360
0361 struct Diagonal2Dense {};
0362
0363 template<> struct AssignmentKind<DenseShape,DiagonalShape> { typedef Diagonal2Dense Kind; };
0364
0365 // Diagonal matrix to Dense assignment
0366 template< typename DstXprType, typename SrcXprType, typename Functor>
0367 struct Assignment<DstXprType, SrcXprType, Functor, Diagonal2Dense>
0368 {
0369   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &/*func*/)
0370   {
0371     Index dstRows = src.rows();
0372     Index dstCols = src.cols();
0373     if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
0374       dst.resize(dstRows, dstCols);
0375
0376     dst.setZero();
0377     dst.diagonal() = src.diagonal();
0378   }
0379
0380   static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &/*func*/)
0381   { dst.diagonal() += src.diagonal(); }
0382
0383   static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &/*func*/)
0384   { dst.diagonal() -= src.diagonal(); }
0385 };
0386
0387 } // namespace internal
0388
0389 } // end namespace Eigen
0390
0391 #endif // EIGEN_DIAGONALMATRIX_H