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 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
derived()39 inline const Derived& derived() const { return *static_cast<const Derived*>(this); }
derived()40 inline Derived& derived() { return *static_cast<Derived*>(this); }
41
toDenseMatrix()42 DenseMatrixType toDenseMatrix() const { return derived(); }
43 template<typename DenseDerived>
44 void evalTo(MatrixBase<DenseDerived> &other) const;
45 template<typename DenseDerived>
addTo(MatrixBase<DenseDerived> & other)46 void addTo(MatrixBase<DenseDerived> &other) const
47 { other.diagonal() += diagonal(); }
48 template<typename DenseDerived>
subTo(MatrixBase<DenseDerived> & other)49 void subTo(MatrixBase<DenseDerived> &other) const
50 { other.diagonal() -= diagonal(); }
51
diagonal()52 inline const DiagonalVectorType& diagonal() const { return derived().diagonal(); }
diagonal()53 inline DiagonalVectorType& diagonal() { return derived().diagonal(); }
54
rows()55 inline Index rows() const { return diagonal().size(); }
cols()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> >
inverse()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
82 template<typename Derived>
83 template<typename DenseDerived>
evalTo(MatrixBase<DenseDerived> & other)84 void 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
104 namespace internal {
105 template<typename _Scalar, int SizeAtCompileTime, int MaxSizeAtCompileTime>
106 struct 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 }
117 template<typename _Scalar, int SizeAtCompileTime, int MaxSizeAtCompileTime>
118 class 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
212 namespace internal {
213 template<typename _DiagonalVectorType>
214 struct 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
230 template<typename _DiagonalVectorType>
231 class 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 **/
259 template<typename Derived>
260 inline const DiagonalWrapper<const Derived>
261 MatrixBase<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 */
274 template<typename Derived>
275 bool 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