-
Notifications
You must be signed in to change notification settings - Fork 0
/
SelectionRule.h
228 lines (195 loc) · 6.28 KB
/
SelectionRule.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
// Copyright (C) 2016-2020 Yixuan Qiu <yixuan.qiu@cos.name>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at https://mozilla.org/MPL/2.0/.
#ifndef SPECTRA_SELECTION_RULE_H
#define SPECTRA_SELECTION_RULE_H
#include <vector> // std::vector
#include <cmath> // std::abs
#include <algorithm> // std::sort
#include <complex> // std::complex
#include <utility> // std::pair
#include <stdexcept> // std::invalid_argument
#include "TypeTraits.h"
namespace Spectra {
///
/// \defgroup Enumerations
///
/// Enumeration types for the selection rule of eigenvalues.
///
///
/// \ingroup Enumerations
///
/// The enumeration of selection rules of desired eigenvalues.
///
enum class SortRule
{
LargestMagn, ///< Select eigenvalues with largest magnitude. Magnitude
///< means the absolute value for real numbers and norm for
///< complex numbers. Applies to both symmetric and general
///< eigen solvers.
LargestReal, ///< Select eigenvalues with largest real part. Only for general eigen solvers.
LargestImag, ///< Select eigenvalues with largest imaginary part (in magnitude). Only for general eigen solvers.
LargestAlge, ///< Select eigenvalues with largest algebraic value, considering
///< any negative sign. Only for symmetric eigen solvers.
SmallestMagn, ///< Select eigenvalues with smallest magnitude. Applies to both symmetric and general
///< eigen solvers.
SmallestReal, ///< Select eigenvalues with smallest real part. Only for general eigen solvers.
SmallestImag, ///< Select eigenvalues with smallest imaginary part (in magnitude). Only for general eigen solvers.
SmallestAlge, ///< Select eigenvalues with smallest algebraic value. Only for symmetric eigen solvers.
BothEnds ///< Select eigenvalues half from each end of the spectrum. When
///< `nev` is odd, compute more from the high end. Only for symmetric eigen solvers.
};
/// \cond
// When comparing eigenvalues, we first calculate the "target" to sort.
// For example, if we want to choose the eigenvalues with
// largest magnitude, the target will be -abs(x).
// The minus sign is due to the fact that std::sort() sorts in ascending order.
// Default target: throw an exception
template <typename Scalar, SortRule Rule>
class SortingTarget
{
public:
static ElemType<Scalar> get(const Scalar& val)
{
using std::abs;
throw std::invalid_argument("incompatible selection rule");
return -abs(val);
}
};
// Specialization for SortRule::LargestMagn
// This covers [float, double, complex] x [SortRule::LargestMagn]
template <typename Scalar>
class SortingTarget<Scalar, SortRule::LargestMagn>
{
public:
static ElemType<Scalar> get(const Scalar& val)
{
using std::abs;
return -abs(val);
}
};
// Specialization for SortRule::LargestReal
// This covers [complex] x [SortRule::LargestReal]
template <typename RealType>
class SortingTarget<std::complex<RealType>, SortRule::LargestReal>
{
public:
static RealType get(const std::complex<RealType>& val)
{
return -val.real();
}
};
// Specialization for SortRule::LargestImag
// This covers [complex] x [SortRule::LargestImag]
template <typename RealType>
class SortingTarget<std::complex<RealType>, SortRule::LargestImag>
{
public:
static RealType get(const std::complex<RealType>& val)
{
using std::abs;
return -abs(val.imag());
}
};
// Specialization for SortRule::LargestAlge
// This covers [float, double] x [SortRule::LargestAlge]
template <typename Scalar>
class SortingTarget<Scalar, SortRule::LargestAlge>
{
public:
static Scalar get(const Scalar& val)
{
return -val;
}
};
// Here SortRule::BothEnds is the same as SortRule::LargestAlge, but
// we need some additional steps, which are done in
// SymEigsSolver.h => retrieve_ritzpair().
// There we move the smallest values to the proper locations.
template <typename Scalar>
class SortingTarget<Scalar, SortRule::BothEnds>
{
public:
static Scalar get(const Scalar& val)
{
return -val;
}
};
// Specialization for SortRule::SmallestMagn
// This covers [float, double, complex] x [SortRule::SmallestMagn]
template <typename Scalar>
class SortingTarget<Scalar, SortRule::SmallestMagn>
{
public:
static ElemType<Scalar> get(const Scalar& val)
{
using std::abs;
return abs(val);
}
};
// Specialization for SortRule::SmallestReal
// This covers [complex] x [SortRule::SmallestReal]
template <typename RealType>
class SortingTarget<std::complex<RealType>, SortRule::SmallestReal>
{
public:
static RealType get(const std::complex<RealType>& val)
{
return val.real();
}
};
// Specialization for SortRule::SmallestImag
// This covers [complex] x [SortRule::SmallestImag]
template <typename RealType>
class SortingTarget<std::complex<RealType>, SortRule::SmallestImag>
{
public:
static RealType get(const std::complex<RealType>& val)
{
using std::abs;
return abs(val.imag());
}
};
// Specialization for SortRule::SmallestAlge
// This covers [float, double] x [SortRule::SmallestAlge]
template <typename Scalar>
class SortingTarget<Scalar, SortRule::SmallestAlge>
{
public:
static Scalar get(const Scalar& val)
{
return val;
}
};
// Sort eigenvalues
template <typename T, SortRule Rule>
class SortEigenvalue
{
private:
using Index = Eigen::Index;
using IndexArray = std::vector<Index>;
const T* m_evals;
IndexArray m_index;
public:
// Sort indices according to the eigenvalues they point to
inline bool operator()(Index i, Index j)
{
return SortingTarget<T, Rule>::get(m_evals[i]) < SortingTarget<T, Rule>::get(m_evals[j]);
}
SortEigenvalue(const T* start, Index size) :
m_evals(start), m_index(size)
{
for (Index i = 0; i < size; i++)
{
m_index[i] = i;
}
std::sort(m_index.begin(), m_index.end(), *this);
}
inline IndexArray index() const { return m_index; }
inline void swap(IndexArray& other) { m_index.swap(other); }
};
/// \endcond
} // namespace Spectra
#endif // SPECTRA_SELECTION_RULE_H