Spectra
GenEigsBase.h
1 // Copyright (C) 2018 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 GEN_EIGS_BASE_H
8 #define GEN_EIGS_BASE_H
9 
10 #include <Eigen/Core>
11 #include <vector> // std::vector
12 #include <cmath> // std::abs, std::pow, std::sqrt
13 #include <algorithm> // std::min, std::copy
14 #include <complex> // std::complex, std::conj, std::norm, std::abs
15 #include <stdexcept> // std::invalid_argument
16 
17 #include "Util/TypeTraits.h"
18 #include "Util/SelectionRule.h"
19 #include "Util/CompInfo.h"
20 #include "Util/SimpleRandom.h"
21 #include "MatOp/internal/ArnoldiOp.h"
22 #include "LinAlg/UpperHessenbergQR.h"
23 #include "LinAlg/DoubleShiftQR.h"
24 #include "LinAlg/UpperHessenbergEigen.h"
25 #include "LinAlg/Arnoldi.h"
26 
27 namespace Spectra {
28 
29 
37 template < typename Scalar,
38  int SelectionRule,
39  typename OpType,
40  typename BOpType >
42 {
43 private:
44  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Matrix;
45  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector;
46  typedef Eigen::Array<Scalar, Eigen::Dynamic, 1> Array;
47  typedef Eigen::Array<bool, Eigen::Dynamic, 1> BoolArray;
48  typedef Eigen::Map<Matrix> MapMat;
49  typedef Eigen::Map<Vector> MapVec;
50  typedef Eigen::Map<const Vector> MapConstVec;
51 
52  typedef std::complex<Scalar> Complex;
53  typedef Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic> ComplexMatrix;
54  typedef Eigen::Matrix<Complex, Eigen::Dynamic, 1> ComplexVector;
55 
56  typedef ArnoldiOp<Scalar, OpType, BOpType> ArnoldiOpType;
57  typedef Arnoldi<Scalar, ArnoldiOpType> ArnoldiFac;
58 
59 protected:
60  OpType* m_op; // object to conduct matrix operation,
61  // e.g. matrix-vector product
62  const int m_n; // dimension of matrix A
63  const int m_nev; // number of eigenvalues requested
64  const int m_ncv; // dimension of Krylov subspace in the Arnoldi method
65  int m_nmatop; // number of matrix operations called
66  int m_niter; // number of restarting iterations
67 
68  ArnoldiFac m_fac; // Arnoldi factorization
69 
70  ComplexVector m_ritz_val; // Ritz values
71  ComplexMatrix m_ritz_vec; // Ritz vectors
72  ComplexVector m_ritz_est; // last row of m_ritz_vec
73 
74 private:
75  BoolArray m_ritz_conv; // indicator of the convergence of Ritz values
76  int m_info; // status of the computation
77 
78  const Scalar m_near_0; // a very small value, but 1.0 / m_near_0 does not overflow
79  // ~= 1e-307 for the "double" type
80  const Scalar m_eps; // the machine precision, ~= 1e-16 for the "double" type
81  const Scalar m_eps23; // m_eps^(2/3), used to test the convergence
82 
83  // Real Ritz values calculated from UpperHessenbergEigen have exact zero imaginary part
84  // Complex Ritz values have exact conjugate pairs
85  // So we use exact tests here
86  static bool is_complex(const Complex& v) { return v.imag() != Scalar(0); }
87  static bool is_conj(const Complex& v1, const Complex& v2) { return v1 == Eigen::numext::conj(v2); }
88 
89  // Implicitly restarted Arnoldi factorization
90  void restart(int k)
91  {
92  using std::norm;
93 
94  if(k >= m_ncv)
95  return;
96 
97  DoubleShiftQR<Scalar> decomp_ds(m_ncv);
98  UpperHessenbergQR<Scalar> decomp_hb(m_ncv);
99  Matrix Q = Matrix::Identity(m_ncv, m_ncv);
100 
101  for(int i = k; i < m_ncv; i++)
102  {
103  if(is_complex(m_ritz_val[i]) && is_conj(m_ritz_val[i], m_ritz_val[i + 1]))
104  {
105  // H - mu * I = Q1 * R1
106  // H <- R1 * Q1 + mu * I = Q1' * H * Q1
107  // H - conj(mu) * I = Q2 * R2
108  // H <- R2 * Q2 + conj(mu) * I = Q2' * H * Q2
109  //
110  // (H - mu * I) * (H - conj(mu) * I) = Q1 * Q2 * R2 * R1 = Q * R
111  const Scalar s = Scalar(2) * m_ritz_val[i].real();
112  const Scalar t = norm(m_ritz_val[i]);
113 
114  decomp_ds.compute(m_fac.matrix_H(), s, t);
115 
116  // Q -> Q * Qi
117  decomp_ds.apply_YQ(Q);
118  // H -> Q'HQ
119  // Matrix Q = Matrix::Identity(m_ncv, m_ncv);
120  // decomp_ds.apply_YQ(Q);
121  // m_fac_H = Q.transpose() * m_fac_H * Q;
122  m_fac.compress_H(decomp_ds);
123 
124  i++;
125  } else {
126  // QR decomposition of H - mu * I, mu is real
127  decomp_hb.compute(m_fac.matrix_H(), m_ritz_val[i].real());
128 
129  // Q -> Q * Qi
130  decomp_hb.apply_YQ(Q);
131  // H -> Q'HQ = RQ + mu * I
132  m_fac.compress_H(decomp_hb);
133  }
134  }
135 
136  m_fac.compress_V(Q);
137  m_fac.factorize_from(k, m_ncv, m_nmatop);
138 
139  retrieve_ritzpair();
140  }
141 
142  // Calculates the number of converged Ritz values
143  int num_converged(Scalar tol)
144  {
145  // thresh = tol * max(m_eps23, abs(theta)), theta for Ritz value
146  Array thresh = tol * m_ritz_val.head(m_nev).array().abs().max(m_eps23);
147  Array resid = m_ritz_est.head(m_nev).array().abs() * m_fac.f_norm();
148  // Converged "wanted" Ritz values
149  m_ritz_conv = (resid < thresh);
150 
151  return m_ritz_conv.cast<int>().sum();
152  }
153 
154  // Returns the adjusted nev for restarting
155  int nev_adjusted(int nconv)
156  {
157  using std::abs;
158 
159  int nev_new = m_nev;
160  for(int i = m_nev; i < m_ncv; i++)
161  if(abs(m_ritz_est[i]) < m_near_0) nev_new++;
162 
163  // Adjust nev_new, according to dnaup2.f line 660~674 in ARPACK
164  nev_new += std::min(nconv, (m_ncv - nev_new) / 2);
165  if(nev_new == 1 && m_ncv >= 6)
166  nev_new = m_ncv / 2;
167  else if(nev_new == 1 && m_ncv > 3)
168  nev_new = 2;
169 
170  if(nev_new > m_ncv - 2)
171  nev_new = m_ncv - 2;
172 
173  // Increase nev by one if ritz_val[nev - 1] and
174  // ritz_val[nev] are conjugate pairs
175  if(is_complex(m_ritz_val[nev_new - 1]) &&
176  is_conj(m_ritz_val[nev_new - 1], m_ritz_val[nev_new]))
177  {
178  nev_new++;
179  }
180 
181  return nev_new;
182  }
183 
184  // Retrieves and sorts Ritz values and Ritz vectors
185  void retrieve_ritzpair()
186  {
187  UpperHessenbergEigen<Scalar> decomp(m_fac.matrix_H());
188  const ComplexVector& evals = decomp.eigenvalues();
189  ComplexMatrix evecs = decomp.eigenvectors();
190 
191  SortEigenvalue<Complex, SelectionRule> sorting(evals.data(), evals.size());
192  std::vector<int> ind = sorting.index();
193 
194  // Copy the Ritz values and vectors to m_ritz_val and m_ritz_vec, respectively
195  for(int i = 0; i < m_ncv; i++)
196  {
197  m_ritz_val[i] = evals[ind[i]];
198  m_ritz_est[i] = evecs(m_ncv - 1, ind[i]);
199  }
200  for(int i = 0; i < m_nev; i++)
201  {
202  m_ritz_vec.col(i).noalias() = evecs.col(ind[i]);
203  }
204  }
205 
206 protected:
207  // Sorts the first nev Ritz pairs in the specified order
208  // This is used to return the final results
209  virtual void sort_ritzpair(int sort_rule)
210  {
211  // First make sure that we have a valid index vector
212  SortEigenvalue<Complex, LARGEST_MAGN> sorting(m_ritz_val.data(), m_nev);
213  std::vector<int> ind = sorting.index();
214 
215  switch(sort_rule)
216  {
217  case LARGEST_MAGN:
218  break;
219  case LARGEST_REAL:
220  {
221  SortEigenvalue<Complex, LARGEST_REAL> sorting(m_ritz_val.data(), m_nev);
222  ind = sorting.index();
223  }
224  break;
225  case LARGEST_IMAG:
226  {
227  SortEigenvalue<Complex, LARGEST_IMAG> sorting(m_ritz_val.data(), m_nev);
228  ind = sorting.index();
229  }
230  break;
231  case SMALLEST_MAGN:
232  {
233  SortEigenvalue<Complex, SMALLEST_MAGN> sorting(m_ritz_val.data(), m_nev);
234  ind = sorting.index();
235  }
236  break;
237  case SMALLEST_REAL:
238  {
239  SortEigenvalue<Complex, SMALLEST_REAL> sorting(m_ritz_val.data(), m_nev);
240  ind = sorting.index();
241  }
242  break;
243  case SMALLEST_IMAG:
244  {
245  SortEigenvalue<Complex, SMALLEST_IMAG> sorting(m_ritz_val.data(), m_nev);
246  ind = sorting.index();
247  }
248  break;
249  default:
250  throw std::invalid_argument("unsupported sorting rule");
251  }
252 
253  ComplexVector new_ritz_val(m_ncv);
254  ComplexMatrix new_ritz_vec(m_ncv, m_nev);
255  BoolArray new_ritz_conv(m_nev);
256 
257  for(int i = 0; i < m_nev; i++)
258  {
259  new_ritz_val[i] = m_ritz_val[ind[i]];
260  new_ritz_vec.col(i).noalias() = m_ritz_vec.col(ind[i]);
261  new_ritz_conv[i] = m_ritz_conv[ind[i]];
262  }
263 
264  m_ritz_val.swap(new_ritz_val);
265  m_ritz_vec.swap(new_ritz_vec);
266  m_ritz_conv.swap(new_ritz_conv);
267  }
268 
269 public:
271 
272  GenEigsBase(OpType* op, BOpType* Bop, int nev, int ncv) :
273  m_op(op),
274  m_n(m_op->rows()),
275  m_nev(nev),
276  m_ncv(ncv > m_n ? m_n : ncv),
277  m_nmatop(0),
278  m_niter(0),
279  m_fac(ArnoldiOpType(op, Bop), m_ncv),
280  m_info(NOT_COMPUTED),
281  m_near_0(TypeTraits<Scalar>::min() * Scalar(10)),
282  m_eps(Eigen::NumTraits<Scalar>::epsilon()),
283  m_eps23(Eigen::numext::pow(m_eps, Scalar(2.0) / 3))
284  {
285  if(nev < 1 || nev > m_n - 2)
286  throw std::invalid_argument("nev must satisfy 1 <= nev <= n - 2, n is the size of matrix");
287 
288  if(ncv < nev + 2 || ncv > m_n)
289  throw std::invalid_argument("ncv must satisfy nev + 2 <= ncv <= n, n is the size of matrix");
290  }
291 
295  virtual ~GenEigsBase() {}
296 
298 
308  void init(const Scalar* init_resid)
309  {
310  // Reset all matrices/vectors to zero
311  m_ritz_val.resize(m_ncv);
312  m_ritz_vec.resize(m_ncv, m_nev);
313  m_ritz_est.resize(m_ncv);
314  m_ritz_conv.resize(m_nev);
315 
316  m_ritz_val.setZero();
317  m_ritz_vec.setZero();
318  m_ritz_est.setZero();
319  m_ritz_conv.setZero();
320 
321  m_nmatop = 0;
322  m_niter = 0;
323 
324  // Initialize the Arnoldi factorization
325  MapConstVec v0(init_resid, m_n);
326  m_fac.init(v0, m_nmatop);
327  }
328 
336  void init()
337  {
338  SimpleRandom<Scalar> rng(0);
339  Vector init_resid = rng.random_vec(m_n);
340  init(init_resid.data());
341  }
342 
363  int compute(int maxit = 1000, Scalar tol = 1e-10, int sort_rule = LARGEST_MAGN)
364  {
365  // The m-step Arnoldi factorization
366  m_fac.factorize_from(1, m_ncv, m_nmatop);
367  retrieve_ritzpair();
368  // Restarting
369  int i, nconv = 0, nev_adj;
370  for(i = 0; i < maxit; i++)
371  {
372  nconv = num_converged(tol);
373  if(nconv >= m_nev)
374  break;
375 
376  nev_adj = nev_adjusted(nconv);
377  restart(nev_adj);
378  }
379  // Sorting results
380  sort_ritzpair(sort_rule);
381 
382  m_niter += i + 1;
383  m_info = (nconv >= m_nev) ? SUCCESSFUL : NOT_CONVERGING;
384 
385  return std::min(m_nev, nconv);
386  }
387 
392  int info() const { return m_info; }
393 
397  int num_iterations() const { return m_niter; }
398 
402  int num_operations() const { return m_nmatop; }
403 
411  ComplexVector eigenvalues() const
412  {
413  const int nconv = m_ritz_conv.cast<int>().sum();
414  ComplexVector res(nconv);
415 
416  if(!nconv)
417  return res;
418 
419  int j = 0;
420  for(int i = 0; i < m_nev; i++)
421  {
422  if(m_ritz_conv[i])
423  {
424  res[j] = m_ritz_val[i];
425  j++;
426  }
427  }
428 
429  return res;
430  }
431 
441  ComplexMatrix eigenvectors(int nvec) const
442  {
443  const int nconv = m_ritz_conv.cast<int>().sum();
444  nvec = std::min(nvec, nconv);
445  ComplexMatrix res(m_n, nvec);
446 
447  if(!nvec)
448  return res;
449 
450  ComplexMatrix ritz_vec_conv(m_ncv, nvec);
451  int j = 0;
452  for(int i = 0; i < m_nev && j < nvec; i++)
453  {
454  if(m_ritz_conv[i])
455  {
456  ritz_vec_conv.col(j).noalias() = m_ritz_vec.col(i);
457  j++;
458  }
459  }
460 
461  res.noalias() = m_fac.matrix_V() * ritz_vec_conv;
462 
463  return res;
464  }
465 
469  ComplexMatrix eigenvectors() const
470  {
471  return eigenvectors(m_nev);
472  }
473 };
474 
475 
476 } // namespace Spectra
477 
478 #endif // GEN_EIGS_BASE_H
int num_iterations() const
Definition: GenEigsBase.h:397
void init(const Scalar *init_resid)
Definition: GenEigsBase.h:308
Computation was successful.
Definition: CompInfo.h:20
Select eigenvalues with smallest imaginary part (in magnitude). Only for general eigen solvers.
Definition: SelectionRule.h:50
ComplexVector eigenvalues() const
Definition: GenEigsBase.h:411
int compute(int maxit=1000, Scalar tol=1e-10, int sort_rule=LARGEST_MAGN)
Definition: GenEigsBase.h:363
ComplexMatrix eigenvectors(int nvec) const
Definition: GenEigsBase.h:441
Select eigenvalues with largest real part. Only for general eigen solvers.
Definition: SelectionRule.h:38
Select eigenvalues with smallest real part. Only for general eigen solvers.
Definition: SelectionRule.h:48
Select eigenvalues with largest imaginary part (in magnitude). Only for general eigen solvers.
Definition: SelectionRule.h:40
int num_operations() const
Definition: GenEigsBase.h:402
ComplexMatrix eigenvectors() const
Definition: GenEigsBase.h:469