Spectra  1.0.1 Header-only C++ Library for Large Scale Eigenvalue Problems
SelectionRule.h
1 // Copyright (C) 2016-2022 Yixuan Qiu <yixuan.qiu@cos.name>
2 //
3 // This Source Code Form is subject to the terms of the Mozilla
4 // Public License v. 2.0. If a copy of the MPL was not distributed
5 // with this file, You can obtain one at https://mozilla.org/MPL/2.0/.
6
7 #ifndef SPECTRA_SELECTION_RULE_H
8 #define SPECTRA_SELECTION_RULE_H
9
10 #include <vector> // std::vector
11 #include <cmath> // std::abs
12 #include <algorithm> // std::sort
13 #include <complex> // std::complex
14 #include <utility> // std::pair
15 #include <stdexcept> // std::invalid_argument
16
17 #include <Eigen/Core>
18 #include "TypeTraits.h"
19
20 namespace Spectra {
21
27
33 enum class SortRule
34 {
35  LargestMagn,
39
40  LargestReal,
41
42  LargestImag,
43
44  LargestAlge,
46
47  SmallestMagn,
49
50  SmallestReal,
51
52  SmallestImag,
53
54  SmallestAlge,
55
56  BothEnds
58 };
59
61
62 // When comparing eigenvalues, we first calculate the "target" to sort.
63 // For example, if we want to choose the eigenvalues with
64 // largest magnitude, the target will be -abs(x).
65 // The minus sign is due to the fact that std::sort() sorts in ascending order.
66
67 // Default target: throw an exception
68 template <typename Scalar, SortRule Rule>
69 class SortingTarget
70 {
71 public:
72  static ElemType<Scalar> get(const Scalar& val)
73  {
74  using std::abs;
75  throw std::invalid_argument("incompatible selection rule");
76  return -abs(val);
77  }
78 };
79
80 // Specialization for SortRule::LargestMagn
81 // This covers [float, double, complex] x [SortRule::LargestMagn]
82 template <typename Scalar>
83 class SortingTarget<Scalar, SortRule::LargestMagn>
84 {
85 public:
86  static ElemType<Scalar> get(const Scalar& val)
87  {
88  using std::abs;
89  return -abs(val);
90  }
91 };
92
93 // Specialization for SortRule::LargestReal
94 // This covers [complex] x [SortRule::LargestReal]
95 template <typename RealType>
96 class SortingTarget<std::complex<RealType>, SortRule::LargestReal>
97 {
98 public:
99  static RealType get(const std::complex<RealType>& val)
100  {
101  return -val.real();
102  }
103 };
104
105 // Specialization for SortRule::LargestImag
106 // This covers [complex] x [SortRule::LargestImag]
107 template <typename RealType>
108 class SortingTarget<std::complex<RealType>, SortRule::LargestImag>
109 {
110 public:
111  static RealType get(const std::complex<RealType>& val)
112  {
113  using std::abs;
114  return -abs(val.imag());
115  }
116 };
117
118 // Specialization for SortRule::LargestAlge
119 // This covers [float, double] x [SortRule::LargestAlge]
120 template <typename Scalar>
121 class SortingTarget<Scalar, SortRule::LargestAlge>
122 {
123 public:
124  static Scalar get(const Scalar& val)
125  {
126  return -val;
127  }
128 };
129
130 // Here SortRule::BothEnds is the same as SortRule::LargestAlge, but
131 // we need some additional steps, which are done in
132 // SymEigsSolver.h => retrieve_ritzpair().
133 // There we move the smallest values to the proper locations.
134 template <typename Scalar>
135 class SortingTarget<Scalar, SortRule::BothEnds>
136 {
137 public:
138  static Scalar get(const Scalar& val)
139  {
140  return -val;
141  }
142 };
143
144 // Specialization for SortRule::SmallestMagn
145 // This covers [float, double, complex] x [SortRule::SmallestMagn]
146 template <typename Scalar>
147 class SortingTarget<Scalar, SortRule::SmallestMagn>
148 {
149 public:
150  static ElemType<Scalar> get(const Scalar& val)
151  {
152  using std::abs;
153  return abs(val);
154  }
155 };
156
157 // Specialization for SortRule::SmallestReal
158 // This covers [complex] x [SortRule::SmallestReal]
159 template <typename RealType>
160 class SortingTarget<std::complex<RealType>, SortRule::SmallestReal>
161 {
162 public:
163  static RealType get(const std::complex<RealType>& val)
164  {
165  return val.real();
166  }
167 };
168
169 // Specialization for SortRule::SmallestImag
170 // This covers [complex] x [SortRule::SmallestImag]
171 template <typename RealType>
172 class SortingTarget<std::complex<RealType>, SortRule::SmallestImag>
173 {
174 public:
175  static RealType get(const std::complex<RealType>& val)
176  {
177  using std::abs;
178  return abs(val.imag());
179  }
180 };
181
182 // Specialization for SortRule::SmallestAlge
183 // This covers [float, double] x [SortRule::SmallestAlge]
184 template <typename Scalar>
185 class SortingTarget<Scalar, SortRule::SmallestAlge>
186 {
187 public:
188  static Scalar get(const Scalar& val)
189  {
190  return val;
191  }
192 };
193
194 // Sort eigenvalues
195 template <typename T, SortRule Rule>
196 class SortEigenvalue
197 {
198 private:
199  using Index = Eigen::Index;
200  using IndexArray = std::vector<Index>;
201
202  const T* m_evals;
203  IndexArray m_index;
204
205 public:
206  // Sort indices according to the eigenvalues they point to
207  inline bool operator()(Index i, Index j)
208  {
209  return SortingTarget<T, Rule>::get(m_evals[i]) < SortingTarget<T, Rule>::get(m_evals[j]);
210  }
211
212  SortEigenvalue(const T* start, Index size) :
213  m_evals(start), m_index(size)
214  {
215  for (Index i = 0; i < size; i++)
216  {
217  m_index[i] = i;
218  }
219  std::sort(m_index.begin(), m_index.end(), *this);
220  }
221
222  inline IndexArray index() const { return m_index; }
223  inline void swap(IndexArray& other) { m_index.swap(other); }
224 };
225
226 // Sort values[:len] according to the selection rule, and return the indices
227 template <typename Scalar>
228 std::vector<Eigen::Index> argsort(SortRule selection, const Eigen::Matrix<Scalar, Eigen::Dynamic, 1>& values, Eigen::Index len)
229 {
230  using Index = Eigen::Index;
231
232  // Sort Ritz values and put the wanted ones at the beginning
233  std::vector<Index> ind;
234  switch (selection)
235  {
237  {
238  SortEigenvalue<Scalar, SortRule::LargestMagn> sorting(values.data(), len);
239  sorting.swap(ind);
240  break;
241  }
242  case SortRule::BothEnds:
244  {
245  SortEigenvalue<Scalar, SortRule::LargestAlge> sorting(values.data(), len);
246  sorting.swap(ind);
247  break;
248  }
250  {
251  SortEigenvalue<Scalar, SortRule::SmallestMagn> sorting(values.data(), len);
252  sorting.swap(ind);
253  break;
254  }
256  {
257  SortEigenvalue<Scalar, SortRule::SmallestAlge> sorting(values.data(), len);
258  sorting.swap(ind);
259  break;
260  }
261  default:
262  throw std::invalid_argument("unsupported selection rule");
263  }
264
265  // For SortRule::BothEnds, the eigenvalues are sorted according to the
266  // SortRule::LargestAlge rule, so we need to move those smallest values to the left
267  // The order would be
268  // Largest => Smallest => 2nd largest => 2nd smallest => ...
269  // We keep this order since the first k values will always be
270  // the wanted collection, no matter k is nev_updated (used in SymEigsBase::restart())
271  // or is nev (used in SymEigsBase::sort_ritzpair())
272  if (selection == SortRule::BothEnds)
273  {
274  std::vector<Index> ind_copy(ind);
275  for (Index i = 0; i < len; i++)
276  {
277  // If i is even, pick values from the left (large values)
278  // If i is odd, pick values from the right (small values)
279  if (i % 2 == 0)
280  ind[i] = ind_copy[i / 2];
281  else
282  ind[i] = ind_copy[len - 1 - i / 2];
283  }
284  }
285
286  return ind;
287 }
288
289 // Default vector length
290 template <typename Scalar>
291 std::vector<Eigen::Index> argsort(SortRule selection, const Eigen::Matrix<Scalar, Eigen::Dynamic, 1>& values)
292 {
293  return argsort<Scalar>(selection, values, values.size());
294 }
295
297
298 } // namespace Spectra
299
300 #endif // SPECTRA_SELECTION_RULE_H
@ 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.
@ SmallestReal
Select eigenvalues with smallest real part. Only for general eigen solvers.
@ LargestImag
Select eigenvalues with largest imaginary part (in magnitude). Only for general eigen solvers.
@ LargestReal
Select eigenvalues with largest real part. Only for general eigen solvers.