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


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:15