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