DiagonalMatrix.h revision c981c48f5bc9aefeffc0bcb0cc3934c2fae179dd 
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 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
82template<typename Derived>
83template<typename DenseDerived>
84void 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
104namespace internal {
105template<typename _Scalar, int SizeAtCompileTime, int MaxSizeAtCompileTime>
106struct 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}
117template<typename _Scalar, int SizeAtCompileTime, int MaxSizeAtCompileTime>
118class 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
212namespace internal {
213template<typename _DiagonalVectorType>
214struct 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
230template<typename _DiagonalVectorType>
231class 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  **/
259template<typename Derived>
260inline const DiagonalWrapper<const Derived>
261MatrixBase<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  */
274template<typename Derived>
275bool 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