Spectra
SelectionRule.h
1 // Copyright (C) 2016-2017 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 http://mozilla.org/MPL/2.0/.
6
7 #ifndef SELECTION_RULE_H
8 #define 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 namespace Spectra {
18
19
25
32 {
34
39
41
43
46
49
51
53
55 };
57
64 {
65  WHICH_LM = 0,
74 };
75
77
78 // Get the element type of a "scalar"
79 // ElemType<double> => double
80 // ElemType< std::complex<double> > => double
81 template <typename T>
82 class ElemType
83 {
84 public:
85  typedef T type;
86 };
87
88 template <typename T>
89 class ElemType< std::complex<T> >
90 {
91 public:
92  typedef T type;
93 };
94
95 // When comparing eigenvalues, we first calculate the "target"
96 // to sort. For example, if we want to choose the eigenvalues with
97 // largest magnitude, the target will be -abs(x).
98 // The minus sign is due to the fact that std::sort() sorts in ascending order.
99
100 // Default target: throw an exception
101 template <typename Scalar, int SelectionRule>
102 class SortingTarget
103 {
104 public:
105  static typename ElemType<Scalar>::type get(const Scalar& val)
106  {
107  using std::abs;
108  throw std::invalid_argument("incompatible selection rule");
109  return -abs(val);
110  }
111 };
112
113 // Specialization for LARGEST_MAGN
114 // This covers [float, double, complex] x [LARGEST_MAGN]
115 template <typename Scalar>
116 class SortingTarget<Scalar, LARGEST_MAGN>
117 {
118 public:
119  static typename ElemType<Scalar>::type get(const Scalar& val)
120  {
121  using std::abs;
122  return -abs(val);
123  }
124 };
125
126 // Specialization for LARGEST_REAL
127 // This covers [complex] x [LARGEST_REAL]
128 template <typename RealType>
129 class SortingTarget<std::complex<RealType>, LARGEST_REAL>
130 {
131 public:
132  static RealType get(const std::complex<RealType>& val)
133  {
134  return -val.real();
135  }
136 };
137
138 // Specialization for LARGEST_IMAG
139 // This covers [complex] x [LARGEST_IMAG]
140 template <typename RealType>
141 class SortingTarget<std::complex<RealType>, LARGEST_IMAG>
142 {
143 public:
144  static RealType get(const std::complex<RealType>& val)
145  {
146  using std::abs;
147  return -abs(val.imag());
148  }
149 };
150
151 // Specialization for LARGEST_ALGE
152 // This covers [float, double] x [LARGEST_ALGE]
153 template <typename Scalar>
154 class SortingTarget<Scalar, LARGEST_ALGE>
155 {
156 public:
157  static Scalar get(const Scalar& val)
158  {
159  return -val;
160  }
161 };
162
163 // Here BOTH_ENDS is the same as LARGEST_ALGE, but
164 // we need some additional steps, which are done in
165 // SymEigsSolver.h => retrieve_ritzpair().
166 // There we move the smallest values to the proper locations.
167 template <typename Scalar>
168 class SortingTarget<Scalar, BOTH_ENDS>
169 {
170 public:
171  static Scalar get(const Scalar& val)
172  {
173  return -val;
174  }
175 };
176
177 // Specialization for SMALLEST_MAGN
178 // This covers [float, double, complex] x [SMALLEST_MAGN]
179 template <typename Scalar>
180 class SortingTarget<Scalar, SMALLEST_MAGN>
181 {
182 public:
183  static typename ElemType<Scalar>::type get(const Scalar& val)
184  {
185  using std::abs;
186  return abs(val);
187  }
188 };
189
190 // Specialization for SMALLEST_REAL
191 // This covers [complex] x [SMALLEST_REAL]
192 template <typename RealType>
193 class SortingTarget<std::complex<RealType>, SMALLEST_REAL>
194 {
195 public:
196  static RealType get(const std::complex<RealType>& val)
197  {
198  return val.real();
199  }
200 };
201
202 // Specialization for SMALLEST_IMAG
203 // This covers [complex] x [SMALLEST_IMAG]
204 template <typename RealType>
205 class SortingTarget<std::complex<RealType>, SMALLEST_IMAG>
206 {
207 public:
208  static RealType get(const std::complex<RealType>& val)
209  {
210  using std::abs;
211  return abs(val.imag());
212  }
213 };
214
215 // Specialization for SMALLEST_ALGE
216 // This covers [float, double] x [SMALLEST_ALGE]
217 template <typename Scalar>
218 class SortingTarget<Scalar, SMALLEST_ALGE>
219 {
220 public:
221  static Scalar get(const Scalar& val)
222  {
223  return val;
224  }
225 };
226
227 // Sort eigenvalues and return the order index
228 template <typename PairType>
229 class PairComparator
230 {
231 public:
232  bool operator() (const PairType& v1, const PairType& v2)
233  {
234  return v1.first < v2.first;
235  }
236 };
237
238 template <typename T, int SelectionRule>
239 class SortEigenvalue
240 {
241 private:
242  typedef typename ElemType<T>::type TargetType; // Type of the sorting target, will be
243  // a floating number type, e.g. "double"
244  typedef std::pair<TargetType, int> PairType; // Type of the sorting pair, including
245  // the sorting target and the index
246
247  std::vector<PairType> pair_sort;
248
249 public:
250  SortEigenvalue(const T* start, int size) :
251  pair_sort(size)
252  {
253  for(int i = 0; i < size; i++)
254  {
255  pair_sort[i].first = SortingTarget<T, SelectionRule>::get(start[i]);
256  pair_sort[i].second = i;
257  }
258  PairComparator<PairType> comp;
259  std::sort(pair_sort.begin(), pair_sort.end(), comp);
260  }
261
262  std::vector<int> index()
263  {
264  std::vector<int> ind(pair_sort.size());
265  for(unsigned int i = 0; i < ind.size(); i++)
266  ind[i] = pair_sort[i].second;
267
268  return ind;
269  }
270 };
271
273
274
275 } // namespace Spectra
276
277 #endif // SELECTION_RULE_H
Alias for SMALLEST_IMAG
Definition: SelectionRule.h:71
Select eigenvalues with smallest imaginary part (in magnitude). Only for general eigen solvers...
Definition: SelectionRule.h:50
Select eigenvalues with smallest algebraic value. Only for symmetric eigen solvers.
Definition: SelectionRule.h:52
Alias for SMALLEST_MAGN
Definition: SelectionRule.h:69
Alias for SMALLEST_REAL
Definition: SelectionRule.h:70
Select eigenvalues with largest real part. Only for general eigen solvers.
Definition: SelectionRule.h:38
Alias for LARGEST_ALGE
Definition: SelectionRule.h:68
Select eigenvalues with smallest real part. Only for general eigen solvers.
Definition: SelectionRule.h:48
SELECT_EIGENVALUE_ALIAS
Definition: SelectionRule.h:63
Select eigenvalues with largest imaginary part (in magnitude). Only for general eigen solvers...
Definition: SelectionRule.h:40
Alias for LARGEST_REAL
Definition: SelectionRule.h:66
Alias for LARGEST_IMAG
Definition: SelectionRule.h:67
Alias for LARGEST_MAGN
Definition: SelectionRule.h:65
Alias for BOTH_ENDS
Definition: SelectionRule.h:73
Alias for SMALLEST_ALGE
Definition: SelectionRule.h:72