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
13namespace 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
31namespace internal {
32template<typename MatrixType, unsigned int UpLo>
33struct traits<SelfAdjointView<MatrixType, UpLo> > : traits<MatrixType>
34{
35 typedef typename ref_selector<MatrixType>::non_const_type MatrixTypeNested;
36 typedef typename remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned;
37 typedef MatrixType ExpressionType;
38 typedef typename MatrixType::PlainObject FullMatrixType;
39 enum {
40 Mode = UpLo | SelfAdjoint,
41 FlagsLvalueBit = is_lvalue<MatrixType>::value ? LvalueBit : 0,
42 Flags = MatrixTypeNestedCleaned::Flags & (HereditaryBits|FlagsLvalueBit)
43 & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit)) // FIXME these flags should be preserved
44 };
45};
46}
47
48
49template<typename _MatrixType, unsigned int UpLo> class SelfAdjointView
50 : public TriangularBase<SelfAdjointView<_MatrixType, UpLo> >
51{
52 public:
53
54 typedef _MatrixType MatrixType;
55 typedef TriangularBase<SelfAdjointView> Base;
56 typedef typename internal::traits<SelfAdjointView>::MatrixTypeNested MatrixTypeNested;
57 typedef typename internal::traits<SelfAdjointView>::MatrixTypeNestedCleaned MatrixTypeNestedCleaned;
58 typedef MatrixTypeNestedCleaned NestedExpression;
59
60 /** \brief The type of coefficients in this matrix */
61 typedef typename internal::traits<SelfAdjointView>::Scalar Scalar;
62 typedef typename MatrixType::StorageIndex StorageIndex;
63 typedef typename internal::remove_all<typename MatrixType::ConjugateReturnType>::type MatrixConjugateReturnType;
64
65 enum {
66 Mode = internal::traits<SelfAdjointView>::Mode,
67 Flags = internal::traits<SelfAdjointView>::Flags,
68 TransposeMode = ((Mode & Upper) ? Lower : 0) | ((Mode & Lower) ? Upper : 0)
69 };
70 typedef typename MatrixType::PlainObject PlainObject;
71
72 EIGEN_DEVICE_FUNC
73 explicit inline SelfAdjointView(MatrixType& matrix) : m_matrix(matrix)
74 {
75 EIGEN_STATIC_ASSERT(UpLo==Lower || UpLo==Upper,SELFADJOINTVIEW_ACCEPTS_UPPER_AND_LOWER_MODE_ONLY);
76 }
77
78 EIGEN_DEVICE_FUNC
79 inline Index rows() const { return m_matrix.rows(); }
80 EIGEN_DEVICE_FUNC
81 inline Index cols() const { return m_matrix.cols(); }
82 EIGEN_DEVICE_FUNC
83 inline Index outerStride() const { return m_matrix.outerStride(); }
84 EIGEN_DEVICE_FUNC
85 inline Index innerStride() const { return m_matrix.innerStride(); }
86
87 /** \sa MatrixBase::coeff()
88 * \warning the coordinates must fit into the referenced triangular part
89 */
90 EIGEN_DEVICE_FUNC
91 inline Scalar coeff(Index row, Index col) const
92 {
93 Base::check_coordinates_internal(row, col);
94 return m_matrix.coeff(row, col);
95 }
96
97 /** \sa MatrixBase::coeffRef()
98 * \warning the coordinates must fit into the referenced triangular part
99 */
100 EIGEN_DEVICE_FUNC
101 inline Scalar& coeffRef(Index row, Index col)
102 {
103 EIGEN_STATIC_ASSERT_LVALUE(SelfAdjointView);
104 Base::check_coordinates_internal(row, col);
105 return m_matrix.coeffRef(row, col);
106 }
107
108 /** \internal */
109 EIGEN_DEVICE_FUNC
110 const MatrixTypeNestedCleaned& _expression() const { return m_matrix; }
111
112 EIGEN_DEVICE_FUNC
113 const MatrixTypeNestedCleaned& nestedExpression() const { return m_matrix; }
114 EIGEN_DEVICE_FUNC
115 MatrixTypeNestedCleaned& nestedExpression() { return m_matrix; }
116
117 /** Efficient triangular matrix times vector/matrix product */
118 template<typename OtherDerived>
119 EIGEN_DEVICE_FUNC
120 const Product<SelfAdjointView,OtherDerived>
121 operator*(const MatrixBase<OtherDerived>& rhs) const
122 {
123 return Product<SelfAdjointView,OtherDerived>(*this, rhs.derived());
124 }
125
126 /** Efficient vector/matrix times triangular matrix product */
127 template<typename OtherDerived> friend
128 EIGEN_DEVICE_FUNC
129 const Product<OtherDerived,SelfAdjointView>
130 operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView& rhs)
131 {
132 return Product<OtherDerived,SelfAdjointView>(lhs.derived(),rhs);
133 }
134
135 friend EIGEN_DEVICE_FUNC
136 const SelfAdjointView<const EIGEN_SCALAR_BINARYOP_EXPR_RETURN_TYPE(Scalar,MatrixType,product),UpLo>
137 operator*(const Scalar& s, const SelfAdjointView& mat)
138 {
139 return (s*mat.nestedExpression()).template selfadjointView<UpLo>();
140 }
141
142 /** Perform a symmetric rank 2 update of the selfadjoint matrix \c *this:
143 * \f$ this = this + \alpha u v^* + conj(\alpha) v u^* \f$
144 * \returns a reference to \c *this
145 *
146 * The vectors \a u and \c v \b must be column vectors, however they can be
147 * a adjoint expression without any overhead. Only the meaningful triangular
148 * part of the matrix is updated, the rest is left unchanged.
149 *
150 * \sa rankUpdate(const MatrixBase<DerivedU>&, Scalar)
151 */
152 template<typename DerivedU, typename DerivedV>
153 EIGEN_DEVICE_FUNC
154 SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, const Scalar& alpha = Scalar(1));
155
156 /** Perform a symmetric rank K update of the selfadjoint matrix \c *this:
157 * \f$ this = this + \alpha ( u u^* ) \f$ where \a u is a vector or matrix.
158 *
159 * \returns a reference to \c *this
160 *
161 * Note that to perform \f$ this = this + \alpha ( u^* u ) \f$ you can simply
162 * call this function with u.adjoint().
163 *
164 * \sa rankUpdate(const MatrixBase<DerivedU>&, const MatrixBase<DerivedV>&, Scalar)
165 */
166 template<typename DerivedU>
167 EIGEN_DEVICE_FUNC
168 SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const Scalar& alpha = Scalar(1));
169
170 /** \returns an expression of a triangular view extracted from the current selfadjoint view of a given triangular part
171 *
172 * The parameter \a TriMode can have the following values: \c #Upper, \c #StrictlyUpper, \c #UnitUpper,
173 * \c #Lower, \c #StrictlyLower, \c #UnitLower.
174 *
175 * If \c TriMode references the same triangular part than \c *this, then this method simply return a \c TriangularView of the nested expression,
176 * otherwise, the nested expression is first transposed, thus returning a \c TriangularView<Transpose<MatrixType>> object.
177 *
178 * \sa MatrixBase::triangularView(), class TriangularView
179 */
180 template<unsigned int TriMode>
181 EIGEN_DEVICE_FUNC
182 typename internal::conditional<(TriMode&(Upper|Lower))==(UpLo&(Upper|Lower)),
183 TriangularView<MatrixType,TriMode>,
184 TriangularView<typename MatrixType::AdjointReturnType,TriMode> >::type
185 triangularView() const
186 {
187 typename internal::conditional<(TriMode&(Upper|Lower))==(UpLo&(Upper|Lower)), MatrixType&, typename MatrixType::ConstTransposeReturnType>::type tmp1(m_matrix);
188 typename internal::conditional<(TriMode&(Upper|Lower))==(UpLo&(Upper|Lower)), MatrixType&, typename MatrixType::AdjointReturnType>::type tmp2(tmp1);
189 return typename internal::conditional<(TriMode&(Upper|Lower))==(UpLo&(Upper|Lower)),
190 TriangularView<MatrixType,TriMode>,
191 TriangularView<typename MatrixType::AdjointReturnType,TriMode> >::type(tmp2);
192 }
193
194 typedef SelfAdjointView<const MatrixConjugateReturnType,UpLo> ConjugateReturnType;
195 /** \sa MatrixBase::conjugate() const */
196 EIGEN_DEVICE_FUNC
197 inline const ConjugateReturnType conjugate() const
198 { return ConjugateReturnType(m_matrix.conjugate()); }
199
200 typedef SelfAdjointView<const typename MatrixType::AdjointReturnType,TransposeMode> AdjointReturnType;
201 /** \sa MatrixBase::adjoint() const */
202 EIGEN_DEVICE_FUNC
203 inline const AdjointReturnType adjoint() const
204 { return AdjointReturnType(m_matrix.adjoint()); }
205
206 typedef SelfAdjointView<typename MatrixType::TransposeReturnType,TransposeMode> TransposeReturnType;
207 /** \sa MatrixBase::transpose() */
208 EIGEN_DEVICE_FUNC
209 inline TransposeReturnType transpose()
210 {
211 EIGEN_STATIC_ASSERT_LVALUE(MatrixType)
212 typename MatrixType::TransposeReturnType tmp(m_matrix);
213 return TransposeReturnType(tmp);
214 }
215
216 typedef SelfAdjointView<const typename MatrixType::ConstTransposeReturnType,TransposeMode> ConstTransposeReturnType;
217 /** \sa MatrixBase::transpose() const */
218 EIGEN_DEVICE_FUNC
219 inline const ConstTransposeReturnType transpose() const
220 {
221 return ConstTransposeReturnType(m_matrix.transpose());
222 }
223
224 /** \returns a const expression of the main diagonal of the matrix \c *this
225 *
226 * This method simply returns the diagonal of the nested expression, thus by-passing the SelfAdjointView decorator.
227 *
228 * \sa MatrixBase::diagonal(), class Diagonal */
229 EIGEN_DEVICE_FUNC
230 typename MatrixType::ConstDiagonalReturnType diagonal() const
231 {
232 return typename MatrixType::ConstDiagonalReturnType(m_matrix);
233 }
234
235/////////// Cholesky module ///////////
236
237 const LLT<PlainObject, UpLo> llt() const;
238 const LDLT<PlainObject, UpLo> ldlt() const;
239
240/////////// Eigenvalue module ///////////
241
242 /** Real part of #Scalar */
243 typedef typename NumTraits<Scalar>::Real RealScalar;
244 /** Return type of eigenvalues() */
245 typedef Matrix<RealScalar, internal::traits<MatrixType>::ColsAtCompileTime, 1> EigenvaluesReturnType;
246
247 EIGEN_DEVICE_FUNC
248 EigenvaluesReturnType eigenvalues() const;
249 EIGEN_DEVICE_FUNC
250 RealScalar operatorNorm() const;
251
252 protected:
253 MatrixTypeNested m_matrix;
254};
255
256
257// template<typename OtherDerived, typename MatrixType, unsigned int UpLo>
258// internal::selfadjoint_matrix_product_returntype<OtherDerived,SelfAdjointView<MatrixType,UpLo> >
259// operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView<MatrixType,UpLo>& rhs)
260// {
261// return internal::matrix_selfadjoint_product_returntype<OtherDerived,SelfAdjointView<MatrixType,UpLo> >(lhs.derived(),rhs);
262// }
263
264// selfadjoint to dense matrix
265
266namespace internal {
267
268// TODO currently a selfadjoint expression has the form SelfAdjointView<.,.>
269// in the future selfadjoint-ness should be defined by the expression traits
270// such that Transpose<SelfAdjointView<.,.> > is valid. (currently TriangularBase::transpose() is overloaded to make it work)
271template<typename MatrixType, unsigned int Mode>
272struct evaluator_traits<SelfAdjointView<MatrixType,Mode> >
273{
274 typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
275 typedef SelfAdjointShape Shape;
276};
277
278template<int UpLo, int SetOpposite, typename DstEvaluatorTypeT, typename SrcEvaluatorTypeT, typename Functor, int Version>
279class triangular_dense_assignment_kernel<UpLo,SelfAdjoint,SetOpposite,DstEvaluatorTypeT,SrcEvaluatorTypeT,Functor,Version>
280 : public generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version>
281{
282protected:
283 typedef generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version> Base;
284 typedef typename Base::DstXprType DstXprType;
285 typedef typename Base::SrcXprType SrcXprType;
286 using Base::m_dst;
287 using Base::m_src;
288 using Base::m_functor;
289public:
290
291 typedef typename Base::DstEvaluatorType DstEvaluatorType;
292 typedef typename Base::SrcEvaluatorType SrcEvaluatorType;
293 typedef typename Base::Scalar Scalar;
294 typedef typename Base::AssignmentTraits AssignmentTraits;
295
296
297 EIGEN_DEVICE_FUNC triangular_dense_assignment_kernel(DstEvaluatorType &dst, const SrcEvaluatorType &src, const Functor &func, DstXprType& dstExpr)
298 : Base(dst, src, func, dstExpr)
299 {}
300
301 EIGEN_DEVICE_FUNC void assignCoeff(Index row, Index col)
302 {
303 eigen_internal_assert(row!=col);
304 Scalar tmp = m_src.coeff(row,col);
305 m_functor.assignCoeff(m_dst.coeffRef(row,col), tmp);
306 m_functor.assignCoeff(m_dst.coeffRef(col,row), numext::conj(tmp));
307 }
308
309 EIGEN_DEVICE_FUNC void assignDiagonalCoeff(Index id)
310 {
311 Base::assignCoeff(id,id);
312 }
313
314 EIGEN_DEVICE_FUNC void assignOppositeCoeff(Index, Index)
315 { eigen_internal_assert(false && "should never be called"); }
316};
317
318} // end namespace internal
319
320/***************************************************************************
321* Implementation of MatrixBase methods
322***************************************************************************/
323
324/** This is the const version of MatrixBase::selfadjointView() */
325template<typename Derived>
326template<unsigned int UpLo>
327typename MatrixBase<Derived>::template ConstSelfAdjointViewReturnType<UpLo>::Type
328MatrixBase<Derived>::selfadjointView() const
329{
330 return typename ConstSelfAdjointViewReturnType<UpLo>::Type(derived());
331}
332
333/** \returns an expression of a symmetric/self-adjoint view extracted from the upper or lower triangular part of the current matrix
334 *
335 * The parameter \a UpLo can be either \c #Upper or \c #Lower
336 *
337 * Example: \include MatrixBase_selfadjointView.cpp
338 * Output: \verbinclude MatrixBase_selfadjointView.out
339 *
340 * \sa class SelfAdjointView
341 */
342template<typename Derived>
343template<unsigned int UpLo>
344typename MatrixBase<Derived>::template SelfAdjointViewReturnType<UpLo>::Type
345MatrixBase<Derived>::selfadjointView()
346{
347 return typename SelfAdjointViewReturnType<UpLo>::Type(derived());
348}
349
350} // end namespace Eigen
351
352#endif // EIGEN_SELFADJOINTMATRIX_H
353