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>::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
88#endif
89
90/** \class DiagonalMatrix
91  * \ingroup Core_Module
92  *
93  * \brief Represents a diagonal matrix with its storage
94  *
95  * \param _Scalar the type of coefficients
96  * \param SizeAtCompileTime the dimension of the matrix, or Dynamic
97  * \param MaxSizeAtCompileTime the dimension of the matrix, or Dynamic. This parameter is optional and defaults
98  *        to SizeAtCompileTime. Most of the time, you do not need to specify it.
99  *
100  * \sa class DiagonalWrapper
101  */
102
103namespace internal {
104template<typename _Scalar, int SizeAtCompileTime, int MaxSizeAtCompileTime>
105struct traits<DiagonalMatrix<_Scalar,SizeAtCompileTime,MaxSizeAtCompileTime> >
106 : traits<Matrix<_Scalar,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
107{
108  typedef Matrix<_Scalar,SizeAtCompileTime,1,0,MaxSizeAtCompileTime,1> DiagonalVectorType;
109  typedef DiagonalShape StorageKind;
110  enum {
111    Flags = LvalueBit | NoPreferredStorageOrderBit
112  };
113};
114}
115template<typename _Scalar, int SizeAtCompileTime, int MaxSizeAtCompileTime>
116class DiagonalMatrix
117  : public DiagonalBase<DiagonalMatrix<_Scalar,SizeAtCompileTime,MaxSizeAtCompileTime> >
118{
119  public:
120    #ifndef EIGEN_PARSED_BY_DOXYGEN
121    typedef typename internal::traits<DiagonalMatrix>::DiagonalVectorType DiagonalVectorType;
122    typedef const DiagonalMatrix& Nested;
123    typedef _Scalar Scalar;
124    typedef typename internal::traits<DiagonalMatrix>::StorageKind StorageKind;
125    typedef typename internal::traits<DiagonalMatrix>::StorageIndex StorageIndex;
126    #endif
127
128  protected:
129
130    DiagonalVectorType m_diagonal;
131
132  public:
133
134    /** const version of diagonal(). */
135    EIGEN_DEVICE_FUNC
136    inline const DiagonalVectorType& diagonal() const { return m_diagonal; }
137    /** \returns a reference to the stored vector of diagonal coefficients. */
138    EIGEN_DEVICE_FUNC
139    inline DiagonalVectorType& diagonal() { return m_diagonal; }
140
141    /** Default constructor without initialization */
142    EIGEN_DEVICE_FUNC
143    inline DiagonalMatrix() {}
144
145    /** Constructs a diagonal matrix with given dimension  */
146    EIGEN_DEVICE_FUNC
147    explicit inline DiagonalMatrix(Index dim) : m_diagonal(dim) {}
148
149    /** 2D constructor. */
150    EIGEN_DEVICE_FUNC
151    inline DiagonalMatrix(const Scalar& x, const Scalar& y) : m_diagonal(x,y) {}
152
153    /** 3D constructor. */
154    EIGEN_DEVICE_FUNC
155    inline DiagonalMatrix(const Scalar& x, const Scalar& y, const Scalar& z) : m_diagonal(x,y,z) {}
156
157    /** Copy constructor. */
158    template<typename OtherDerived>
159    EIGEN_DEVICE_FUNC
160    inline DiagonalMatrix(const DiagonalBase<OtherDerived>& other) : m_diagonal(other.diagonal()) {}
161
162    #ifndef EIGEN_PARSED_BY_DOXYGEN
163    /** copy constructor. prevent a default copy constructor from hiding the other templated constructor */
164    inline DiagonalMatrix(const DiagonalMatrix& other) : m_diagonal(other.diagonal()) {}
165    #endif
166
167    /** generic constructor from expression of the diagonal coefficients */
168    template<typename OtherDerived>
169    EIGEN_DEVICE_FUNC
170    explicit inline DiagonalMatrix(const MatrixBase<OtherDerived>& other) : m_diagonal(other)
171    {}
172
173    /** Copy operator. */
174    template<typename OtherDerived>
175    EIGEN_DEVICE_FUNC
176    DiagonalMatrix& operator=(const DiagonalBase<OtherDerived>& other)
177    {
178      m_diagonal = other.diagonal();
179      return *this;
180    }
181
182    #ifndef EIGEN_PARSED_BY_DOXYGEN
183    /** This is a special case of the templated operator=. Its purpose is to
184      * prevent a default operator= from hiding the templated operator=.
185      */
186    EIGEN_DEVICE_FUNC
187    DiagonalMatrix& operator=(const DiagonalMatrix& other)
188    {
189      m_diagonal = other.diagonal();
190      return *this;
191    }
192    #endif
193
194    /** Resizes to given size. */
195    EIGEN_DEVICE_FUNC
196    inline void resize(Index size) { m_diagonal.resize(size); }
197    /** Sets all coefficients to zero. */
198    EIGEN_DEVICE_FUNC
199    inline void setZero() { m_diagonal.setZero(); }
200    /** Resizes and sets all coefficients to zero. */
201    EIGEN_DEVICE_FUNC
202    inline void setZero(Index size) { m_diagonal.setZero(size); }
203    /** Sets this matrix to be the identity matrix of the current size. */
204    EIGEN_DEVICE_FUNC
205    inline void setIdentity() { m_diagonal.setOnes(); }
206    /** Sets this matrix to be the identity matrix of the given size. */
207    EIGEN_DEVICE_FUNC
208    inline void setIdentity(Index size) { m_diagonal.setOnes(size); }
209};
210
211/** \class DiagonalWrapper
212  * \ingroup Core_Module
213  *
214  * \brief Expression of a diagonal matrix
215  *
216  * \param _DiagonalVectorType the type of the vector of diagonal coefficients
217  *
218  * This class is an expression of a diagonal matrix, but not storing its own vector of diagonal coefficients,
219  * instead wrapping an existing vector expression. It is the return type of MatrixBase::asDiagonal()
220  * and most of the time this is the only way that it is used.
221  *
222  * \sa class DiagonalMatrix, class DiagonalBase, MatrixBase::asDiagonal()
223  */
224
225namespace internal {
226template<typename _DiagonalVectorType>
227struct traits<DiagonalWrapper<_DiagonalVectorType> >
228{
229  typedef _DiagonalVectorType DiagonalVectorType;
230  typedef typename DiagonalVectorType::Scalar Scalar;
231  typedef typename DiagonalVectorType::StorageIndex StorageIndex;
232  typedef DiagonalShape StorageKind;
233  typedef typename traits<DiagonalVectorType>::XprKind XprKind;
234  enum {
235    RowsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
236    ColsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
237    MaxRowsAtCompileTime = DiagonalVectorType::MaxSizeAtCompileTime,
238    MaxColsAtCompileTime = DiagonalVectorType::MaxSizeAtCompileTime,
239    Flags =  (traits<DiagonalVectorType>::Flags & LvalueBit) | NoPreferredStorageOrderBit
240  };
241};
242}
243
244template<typename _DiagonalVectorType>
245class DiagonalWrapper
246  : public DiagonalBase<DiagonalWrapper<_DiagonalVectorType> >, internal::no_assignment_operator
247{
248  public:
249    #ifndef EIGEN_PARSED_BY_DOXYGEN
250    typedef _DiagonalVectorType DiagonalVectorType;
251    typedef DiagonalWrapper Nested;
252    #endif
253
254    /** Constructor from expression of diagonal coefficients to wrap. */
255    EIGEN_DEVICE_FUNC
256    explicit inline DiagonalWrapper(DiagonalVectorType& a_diagonal) : m_diagonal(a_diagonal) {}
257
258    /** \returns a const reference to the wrapped expression of diagonal coefficients. */
259    EIGEN_DEVICE_FUNC
260    const DiagonalVectorType& diagonal() const { return m_diagonal; }
261
262  protected:
263    typename DiagonalVectorType::Nested m_diagonal;
264};
265
266/** \returns a pseudo-expression of a diagonal matrix with *this as vector of diagonal coefficients
267  *
268  * \only_for_vectors
269  *
270  * Example: \include MatrixBase_asDiagonal.cpp
271  * Output: \verbinclude MatrixBase_asDiagonal.out
272  *
273  * \sa class DiagonalWrapper, class DiagonalMatrix, diagonal(), isDiagonal()
274  **/
275template<typename Derived>
276inline const DiagonalWrapper<const Derived>
277MatrixBase<Derived>::asDiagonal() const
278{
279  return DiagonalWrapper<const Derived>(derived());
280}
281
282/** \returns true if *this is approximately equal to a diagonal matrix,
283  *          within the precision given by \a prec.
284  *
285  * Example: \include MatrixBase_isDiagonal.cpp
286  * Output: \verbinclude MatrixBase_isDiagonal.out
287  *
288  * \sa asDiagonal()
289  */
290template<typename Derived>
291bool MatrixBase<Derived>::isDiagonal(const RealScalar& prec) const
292{
293  if(cols() != rows()) return false;
294  RealScalar maxAbsOnDiagonal = static_cast<RealScalar>(-1);
295  for(Index j = 0; j < cols(); ++j)
296  {
297    RealScalar absOnDiagonal = numext::abs(coeff(j,j));
298    if(absOnDiagonal > maxAbsOnDiagonal) maxAbsOnDiagonal = absOnDiagonal;
299  }
300  for(Index j = 0; j < cols(); ++j)
301    for(Index i = 0; i < j; ++i)
302    {
303      if(!internal::isMuchSmallerThan(coeff(i, j), maxAbsOnDiagonal, prec)) return false;
304      if(!internal::isMuchSmallerThan(coeff(j, i), maxAbsOnDiagonal, prec)) return false;
305    }
306  return true;
307}
308
309namespace internal {
310
311template<> struct storage_kind_to_shape<DiagonalShape> { typedef DiagonalShape Shape; };
312
313struct Diagonal2Dense {};
314
315template<> struct AssignmentKind<DenseShape,DiagonalShape> { typedef Diagonal2Dense Kind; };
316
317// Diagonal matrix to Dense assignment
318template< typename DstXprType, typename SrcXprType, typename Functor>
319struct Assignment<DstXprType, SrcXprType, Functor, Diagonal2Dense>
320{
321  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &/*func*/)
322  {
323    Index dstRows = src.rows();
324    Index dstCols = src.cols();
325    if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
326      dst.resize(dstRows, dstCols);
327
328    dst.setZero();
329    dst.diagonal() = src.diagonal();
330  }
331
332  static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &/*func*/)
333  { dst.diagonal() += src.diagonal(); }
334
335  static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &/*func*/)
336  { dst.diagonal() -= src.diagonal(); }
337};
338
339} // namespace internal
340
341} // end namespace Eigen
342
343#endif // EIGEN_DIAGONALMATRIX_H
344