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 |
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 | |
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 | |