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
14namespace Eigen {
15
16#ifndef EIGEN_PARSED_BY_DOXYGEN
17template<typename Derived>
18class 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
99template<typename Derived>
100template<typename DenseDerived>
101void 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
121namespace internal {
122template<typename _Scalar, int SizeAtCompileTime, int MaxSizeAtCompileTime>
123struct 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}
134template<typename _Scalar, int SizeAtCompileTime, int MaxSizeAtCompileTime>
135class 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
229namespace internal {
230template<typename _DiagonalVectorType>
231struct 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
247template<typename _DiagonalVectorType>
248class 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  **/
276template<typename Derived>
277inline const DiagonalWrapper<const Derived>
278MatrixBase<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  */
291template<typename Derived>
292bool 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