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 
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 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
derived()41     inline const Derived& derived() const { return *static_cast<const Derived*>(this); }
42     EIGEN_DEVICE_FUNC
derived()43     inline Derived& derived() { return *static_cast<Derived*>(this); }
44 
45     EIGEN_DEVICE_FUNC
toDenseMatrix()46     DenseMatrixType toDenseMatrix() const { return derived(); }
47 
48     EIGEN_DEVICE_FUNC
diagonal()49     inline const DiagonalVectorType& diagonal() const { return derived().diagonal(); }
50     EIGEN_DEVICE_FUNC
diagonal()51     inline DiagonalVectorType& diagonal() { return derived().diagonal(); }
52 
53     EIGEN_DEVICE_FUNC
rows()54     inline Index rows() const { return diagonal().size(); }
55     EIGEN_DEVICE_FUNC
cols()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
inverse()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 
103 namespace internal {
104 template<typename _Scalar, int SizeAtCompileTime, int MaxSizeAtCompileTime>
105 struct 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 }
115 template<typename _Scalar, int SizeAtCompileTime, int MaxSizeAtCompileTime>
116 class 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 
225 namespace internal {
226 template<typename _DiagonalVectorType>
227 struct 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 
244 template<typename _DiagonalVectorType>
245 class 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   **/
275 template<typename Derived>
276 inline const DiagonalWrapper<const Derived>
277 MatrixBase<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   */
290 template<typename Derived>
291 bool 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 
309 namespace internal {
310 
311 template<> struct storage_kind_to_shape<DiagonalShape> { typedef DiagonalShape Shape; };
312 
313 struct Diagonal2Dense {};
314 
315 template<> struct AssignmentKind<DenseShape,DiagonalShape> { typedef Diagonal2Dense Kind; };
316 
317 // Diagonal matrix to Dense assignment
318 template< typename DstXprType, typename SrcXprType, typename Functor>
319 struct 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