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>::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
derived()40 inline const Derived& derived() const { return *static_cast<const Derived*>(this); }
derived()41 inline Derived& derived() { return *static_cast<Derived*>(this); }
42
toDenseMatrix()43 DenseMatrixType toDenseMatrix() const { return derived(); }
44 template<typename DenseDerived>
45 void evalTo(MatrixBase<DenseDerived> &other) const;
46 template<typename DenseDerived>
addTo(MatrixBase<DenseDerived> & other)47 void addTo(MatrixBase<DenseDerived> &other) const
48 { other.diagonal() += diagonal(); }
49 template<typename DenseDerived>
subTo(MatrixBase<DenseDerived> & other)50 void subTo(MatrixBase<DenseDerived> &other) const
51 { other.diagonal() -= diagonal(); }
52
diagonal()53 inline const DiagonalVectorType& diagonal() const { return derived().diagonal(); }
diagonal()54 inline DiagonalVectorType& diagonal() { return derived().diagonal(); }
55
rows()56 inline Index rows() const { return diagonal().size(); }
cols()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> >
inverse()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
99 template<typename Derived>
100 template<typename DenseDerived>
evalTo(MatrixBase<DenseDerived> & other)101 void 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
121 namespace internal {
122 template<typename _Scalar, int SizeAtCompileTime, int MaxSizeAtCompileTime>
123 struct 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 }
134 template<typename _Scalar, int SizeAtCompileTime, int MaxSizeAtCompileTime>
135 class 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
229 namespace internal {
230 template<typename _DiagonalVectorType>
231 struct 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
247 template<typename _DiagonalVectorType>
248 class 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 **/
276 template<typename Derived>
277 inline const DiagonalWrapper<const Derived>
278 MatrixBase<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 */
291 template<typename Derived>
292 bool 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