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 | //
|
---|
6 | // This Source Code Form is subject to the terms of the Mozilla
|
---|
7 | // Public License v. 2.0. If a copy of the MPL was not distributed
|
---|
8 | // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
---|
9 |
|
---|
10 | #ifndef EIGEN_SELFADJOINTMATRIX_H
|
---|
11 | #define EIGEN_SELFADJOINTMATRIX_H
|
---|
12 |
|
---|
13 | namespace Eigen {
|
---|
14 |
|
---|
15 | /** \class SelfAdjointView
|
---|
16 | * \ingroup Core_Module
|
---|
17 | *
|
---|
18 | *
|
---|
19 | * \brief Expression of a selfadjoint matrix from a triangular part of a dense matrix
|
---|
20 | *
|
---|
21 | * \param MatrixType the type of the dense matrix storing the coefficients
|
---|
22 | * \param TriangularPart can be either \c #Lower or \c #Upper
|
---|
23 | *
|
---|
24 | * This class is an expression of a sefladjoint matrix from a triangular part of a matrix
|
---|
25 | * with given dense storage of the coefficients. It is the return type of MatrixBase::selfadjointView()
|
---|
26 | * and most of the time this is the only way that it is used.
|
---|
27 | *
|
---|
28 | * \sa class TriangularBase, MatrixBase::selfadjointView()
|
---|
29 | */
|
---|
30 |
|
---|
31 | namespace internal {
|
---|
32 | template<typename MatrixType, unsigned int UpLo>
|
---|
33 | struct traits<SelfAdjointView<MatrixType, UpLo> > : traits<MatrixType>
|
---|
34 | {
|
---|
35 | typedef typename nested<MatrixType>::type MatrixTypeNested;
|
---|
36 | typedef typename remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned;
|
---|
37 | typedef MatrixType ExpressionType;
|
---|
38 | typedef typename MatrixType::PlainObject DenseMatrixType;
|
---|
39 | enum {
|
---|
40 | Mode = UpLo | SelfAdjoint,
|
---|
41 | Flags = MatrixTypeNestedCleaned::Flags & (HereditaryBits)
|
---|
42 | & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit)), // FIXME these flags should be preserved
|
---|
43 | CoeffReadCost = MatrixTypeNestedCleaned::CoeffReadCost
|
---|
44 | };
|
---|
45 | };
|
---|
46 | }
|
---|
47 |
|
---|
48 | template <typename Lhs, int LhsMode, bool LhsIsVector,
|
---|
49 | typename Rhs, int RhsMode, bool RhsIsVector>
|
---|
50 | struct SelfadjointProductMatrix;
|
---|
51 |
|
---|
52 | // FIXME could also be called SelfAdjointWrapper to be consistent with DiagonalWrapper ??
|
---|
53 | template<typename MatrixType, unsigned int UpLo> class SelfAdjointView
|
---|
54 | : public TriangularBase<SelfAdjointView<MatrixType, UpLo> >
|
---|
55 | {
|
---|
56 | public:
|
---|
57 |
|
---|
58 | typedef TriangularBase<SelfAdjointView> Base;
|
---|
59 | typedef typename internal::traits<SelfAdjointView>::MatrixTypeNested MatrixTypeNested;
|
---|
60 | typedef typename internal::traits<SelfAdjointView>::MatrixTypeNestedCleaned MatrixTypeNestedCleaned;
|
---|
61 |
|
---|
62 | /** \brief The type of coefficients in this matrix */
|
---|
63 | typedef typename internal::traits<SelfAdjointView>::Scalar Scalar;
|
---|
64 |
|
---|
65 | typedef typename MatrixType::Index Index;
|
---|
66 |
|
---|
67 | enum {
|
---|
68 | Mode = internal::traits<SelfAdjointView>::Mode
|
---|
69 | };
|
---|
70 | typedef typename MatrixType::PlainObject PlainObject;
|
---|
71 |
|
---|
72 | inline SelfAdjointView(MatrixType& matrix) : m_matrix(matrix)
|
---|
73 | {}
|
---|
74 |
|
---|
75 | inline Index rows() const { return m_matrix.rows(); }
|
---|
76 | inline Index cols() const { return m_matrix.cols(); }
|
---|
77 | inline Index outerStride() const { return m_matrix.outerStride(); }
|
---|
78 | inline Index innerStride() const { return m_matrix.innerStride(); }
|
---|
79 |
|
---|
80 | /** \sa MatrixBase::coeff()
|
---|
81 | * \warning the coordinates must fit into the referenced triangular part
|
---|
82 | */
|
---|
83 | inline Scalar coeff(Index row, Index col) const
|
---|
84 | {
|
---|
85 | Base::check_coordinates_internal(row, col);
|
---|
86 | return m_matrix.coeff(row, col);
|
---|
87 | }
|
---|
88 |
|
---|
89 | /** \sa MatrixBase::coeffRef()
|
---|
90 | * \warning the coordinates must fit into the referenced triangular part
|
---|
91 | */
|
---|
92 | inline Scalar& coeffRef(Index row, Index col)
|
---|
93 | {
|
---|
94 | Base::check_coordinates_internal(row, col);
|
---|
95 | return m_matrix.const_cast_derived().coeffRef(row, col);
|
---|
96 | }
|
---|
97 |
|
---|
98 | /** \internal */
|
---|
99 | const MatrixTypeNestedCleaned& _expression() const { return m_matrix; }
|
---|
100 |
|
---|
101 | const MatrixTypeNestedCleaned& nestedExpression() const { return m_matrix; }
|
---|
102 | MatrixTypeNestedCleaned& nestedExpression() { return *const_cast<MatrixTypeNestedCleaned*>(&m_matrix); }
|
---|
103 |
|
---|
104 | /** Efficient self-adjoint matrix times vector/matrix product */
|
---|
105 | template<typename OtherDerived>
|
---|
106 | SelfadjointProductMatrix<MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime>
|
---|
107 | operator*(const MatrixBase<OtherDerived>& rhs) const
|
---|
108 | {
|
---|
109 | return SelfadjointProductMatrix
|
---|
110 | <MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime>
|
---|
111 | (m_matrix, rhs.derived());
|
---|
112 | }
|
---|
113 |
|
---|
114 | /** Efficient vector/matrix times self-adjoint matrix product */
|
---|
115 | template<typename OtherDerived> friend
|
---|
116 | SelfadjointProductMatrix<OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false>
|
---|
117 | operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView& rhs)
|
---|
118 | {
|
---|
119 | return SelfadjointProductMatrix
|
---|
120 | <OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false>
|
---|
121 | (lhs.derived(),rhs.m_matrix);
|
---|
122 | }
|
---|
123 |
|
---|
124 | /** Perform a symmetric rank 2 update of the selfadjoint matrix \c *this:
|
---|
125 | * \f$ this = this + \alpha u v^* + conj(\alpha) v u^* \f$
|
---|
126 | * \returns a reference to \c *this
|
---|
127 | *
|
---|
128 | * The vectors \a u and \c v \b must be column vectors, however they can be
|
---|
129 | * a adjoint expression without any overhead. Only the meaningful triangular
|
---|
130 | * part of the matrix is updated, the rest is left unchanged.
|
---|
131 | *
|
---|
132 | * \sa rankUpdate(const MatrixBase<DerivedU>&, Scalar)
|
---|
133 | */
|
---|
134 | template<typename DerivedU, typename DerivedV>
|
---|
135 | SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, const Scalar& alpha = Scalar(1));
|
---|
136 |
|
---|
137 | /** Perform a symmetric rank K update of the selfadjoint matrix \c *this:
|
---|
138 | * \f$ this = this + \alpha ( u u^* ) \f$ where \a u is a vector or matrix.
|
---|
139 | *
|
---|
140 | * \returns a reference to \c *this
|
---|
141 | *
|
---|
142 | * Note that to perform \f$ this = this + \alpha ( u^* u ) \f$ you can simply
|
---|
143 | * call this function with u.adjoint().
|
---|
144 | *
|
---|
145 | * \sa rankUpdate(const MatrixBase<DerivedU>&, const MatrixBase<DerivedV>&, Scalar)
|
---|
146 | */
|
---|
147 | template<typename DerivedU>
|
---|
148 | SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const Scalar& alpha = Scalar(1));
|
---|
149 |
|
---|
150 | /////////// Cholesky module ///////////
|
---|
151 |
|
---|
152 | const LLT<PlainObject, UpLo> llt() const;
|
---|
153 | const LDLT<PlainObject, UpLo> ldlt() const;
|
---|
154 |
|
---|
155 | /////////// Eigenvalue module ///////////
|
---|
156 |
|
---|
157 | /** Real part of #Scalar */
|
---|
158 | typedef typename NumTraits<Scalar>::Real RealScalar;
|
---|
159 | /** Return type of eigenvalues() */
|
---|
160 | typedef Matrix<RealScalar, internal::traits<MatrixType>::ColsAtCompileTime, 1> EigenvaluesReturnType;
|
---|
161 |
|
---|
162 | EigenvaluesReturnType eigenvalues() const;
|
---|
163 | RealScalar operatorNorm() const;
|
---|
164 |
|
---|
165 | #ifdef EIGEN2_SUPPORT
|
---|
166 | template<typename OtherDerived>
|
---|
167 | SelfAdjointView& operator=(const MatrixBase<OtherDerived>& other)
|
---|
168 | {
|
---|
169 | enum {
|
---|
170 | OtherPart = UpLo == Upper ? StrictlyLower : StrictlyUpper
|
---|
171 | };
|
---|
172 | m_matrix.const_cast_derived().template triangularView<UpLo>() = other;
|
---|
173 | m_matrix.const_cast_derived().template triangularView<OtherPart>() = other.adjoint();
|
---|
174 | return *this;
|
---|
175 | }
|
---|
176 | template<typename OtherMatrixType, unsigned int OtherMode>
|
---|
177 | SelfAdjointView& operator=(const TriangularView<OtherMatrixType, OtherMode>& other)
|
---|
178 | {
|
---|
179 | enum {
|
---|
180 | OtherPart = UpLo == Upper ? StrictlyLower : StrictlyUpper
|
---|
181 | };
|
---|
182 | m_matrix.const_cast_derived().template triangularView<UpLo>() = other.toDenseMatrix();
|
---|
183 | m_matrix.const_cast_derived().template triangularView<OtherPart>() = other.toDenseMatrix().adjoint();
|
---|
184 | return *this;
|
---|
185 | }
|
---|
186 | #endif
|
---|
187 |
|
---|
188 | protected:
|
---|
189 | MatrixTypeNested m_matrix;
|
---|
190 | };
|
---|
191 |
|
---|
192 |
|
---|
193 | // template<typename OtherDerived, typename MatrixType, unsigned int UpLo>
|
---|
194 | // internal::selfadjoint_matrix_product_returntype<OtherDerived,SelfAdjointView<MatrixType,UpLo> >
|
---|
195 | // operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView<MatrixType,UpLo>& rhs)
|
---|
196 | // {
|
---|
197 | // return internal::matrix_selfadjoint_product_returntype<OtherDerived,SelfAdjointView<MatrixType,UpLo> >(lhs.derived(),rhs);
|
---|
198 | // }
|
---|
199 |
|
---|
200 | // selfadjoint to dense matrix
|
---|
201 |
|
---|
202 | namespace internal {
|
---|
203 |
|
---|
204 | template<typename Derived1, typename Derived2, int UnrollCount, bool ClearOpposite>
|
---|
205 | struct triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Upper), UnrollCount, ClearOpposite>
|
---|
206 | {
|
---|
207 | enum {
|
---|
208 | col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
|
---|
209 | row = (UnrollCount-1) % Derived1::RowsAtCompileTime
|
---|
210 | };
|
---|
211 |
|
---|
212 | static inline void run(Derived1 &dst, const Derived2 &src)
|
---|
213 | {
|
---|
214 | triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Upper), UnrollCount-1, ClearOpposite>::run(dst, src);
|
---|
215 |
|
---|
216 | if(row == col)
|
---|
217 | dst.coeffRef(row, col) = numext::real(src.coeff(row, col));
|
---|
218 | else if(row < col)
|
---|
219 | dst.coeffRef(col, row) = numext::conj(dst.coeffRef(row, col) = src.coeff(row, col));
|
---|
220 | }
|
---|
221 | };
|
---|
222 |
|
---|
223 | template<typename Derived1, typename Derived2, bool ClearOpposite>
|
---|
224 | struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Upper, 0, ClearOpposite>
|
---|
225 | {
|
---|
226 | static inline void run(Derived1 &, const Derived2 &) {}
|
---|
227 | };
|
---|
228 |
|
---|
229 | template<typename Derived1, typename Derived2, int UnrollCount, bool ClearOpposite>
|
---|
230 | struct triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Lower), UnrollCount, ClearOpposite>
|
---|
231 | {
|
---|
232 | enum {
|
---|
233 | col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
|
---|
234 | row = (UnrollCount-1) % Derived1::RowsAtCompileTime
|
---|
235 | };
|
---|
236 |
|
---|
237 | static inline void run(Derived1 &dst, const Derived2 &src)
|
---|
238 | {
|
---|
239 | triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Lower), UnrollCount-1, ClearOpposite>::run(dst, src);
|
---|
240 |
|
---|
241 | if(row == col)
|
---|
242 | dst.coeffRef(row, col) = numext::real(src.coeff(row, col));
|
---|
243 | else if(row > col)
|
---|
244 | dst.coeffRef(col, row) = numext::conj(dst.coeffRef(row, col) = src.coeff(row, col));
|
---|
245 | }
|
---|
246 | };
|
---|
247 |
|
---|
248 | template<typename Derived1, typename Derived2, bool ClearOpposite>
|
---|
249 | struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Lower, 0, ClearOpposite>
|
---|
250 | {
|
---|
251 | static inline void run(Derived1 &, const Derived2 &) {}
|
---|
252 | };
|
---|
253 |
|
---|
254 | template<typename Derived1, typename Derived2, bool ClearOpposite>
|
---|
255 | struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Upper, Dynamic, ClearOpposite>
|
---|
256 | {
|
---|
257 | typedef typename Derived1::Index Index;
|
---|
258 | static inline void run(Derived1 &dst, const Derived2 &src)
|
---|
259 | {
|
---|
260 | for(Index j = 0; j < dst.cols(); ++j)
|
---|
261 | {
|
---|
262 | for(Index i = 0; i < j; ++i)
|
---|
263 | {
|
---|
264 | dst.copyCoeff(i, j, src);
|
---|
265 | dst.coeffRef(j,i) = numext::conj(dst.coeff(i,j));
|
---|
266 | }
|
---|
267 | dst.copyCoeff(j, j, src);
|
---|
268 | }
|
---|
269 | }
|
---|
270 | };
|
---|
271 |
|
---|
272 | template<typename Derived1, typename Derived2, bool ClearOpposite>
|
---|
273 | struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Lower, Dynamic, ClearOpposite>
|
---|
274 | {
|
---|
275 | static inline void run(Derived1 &dst, const Derived2 &src)
|
---|
276 | {
|
---|
277 | typedef typename Derived1::Index Index;
|
---|
278 | for(Index i = 0; i < dst.rows(); ++i)
|
---|
279 | {
|
---|
280 | for(Index j = 0; j < i; ++j)
|
---|
281 | {
|
---|
282 | dst.copyCoeff(i, j, src);
|
---|
283 | dst.coeffRef(j,i) = numext::conj(dst.coeff(i,j));
|
---|
284 | }
|
---|
285 | dst.copyCoeff(i, i, src);
|
---|
286 | }
|
---|
287 | }
|
---|
288 | };
|
---|
289 |
|
---|
290 | } // end namespace internal
|
---|
291 |
|
---|
292 | /***************************************************************************
|
---|
293 | * Implementation of MatrixBase methods
|
---|
294 | ***************************************************************************/
|
---|
295 |
|
---|
296 | template<typename Derived>
|
---|
297 | template<unsigned int UpLo>
|
---|
298 | typename MatrixBase<Derived>::template ConstSelfAdjointViewReturnType<UpLo>::Type
|
---|
299 | MatrixBase<Derived>::selfadjointView() const
|
---|
300 | {
|
---|
301 | return derived();
|
---|
302 | }
|
---|
303 |
|
---|
304 | template<typename Derived>
|
---|
305 | template<unsigned int UpLo>
|
---|
306 | typename MatrixBase<Derived>::template SelfAdjointViewReturnType<UpLo>::Type
|
---|
307 | MatrixBase<Derived>::selfadjointView()
|
---|
308 | {
|
---|
309 | return derived();
|
---|
310 | }
|
---|
311 |
|
---|
312 | } // end namespace Eigen
|
---|
313 |
|
---|
314 | #endif // EIGEN_SELFADJOINTMATRIX_H
|
---|