Home | History | Annotate | Download | only in Core
      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 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 
     39     inline const Derived& derived() const { return *static_cast<const Derived*>(this); }
     40     inline Derived& derived() { return *static_cast<Derived*>(this); }
     41 
     42     DenseMatrixType toDenseMatrix() const { return derived(); }
     43     template<typename DenseDerived>
     44     void evalTo(MatrixBase<DenseDerived> &other) const;
     45     template<typename DenseDerived>
     46     void addTo(MatrixBase<DenseDerived> &other) const
     47     { other.diagonal() += diagonal(); }
     48     template<typename DenseDerived>
     49     void subTo(MatrixBase<DenseDerived> &other) const
     50     { other.diagonal() -= diagonal(); }
     51 
     52     inline const DiagonalVectorType& diagonal() const { return derived().diagonal(); }
     53     inline DiagonalVectorType& diagonal() { return derived().diagonal(); }
     54 
     55     inline Index rows() const { return diagonal().size(); }
     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> >
     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>
     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