Spectra  1.0.0
Header-only C++ Library for Large Scale Eigenvalue Problems
JDSymEigsBase.h
1 // Copyright (C) 2020 Netherlands eScience Center <J.Wehner@esciencecenter.nl>
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_JD_SYM_EIGS_BASE_H
8 #define SPECTRA_JD_SYM_EIGS_BASE_H
9 
10 #include <Eigen/Core>
11 #include <vector> // std::vector
12 #include <cmath> // std::abs, std::pow
13 #include <algorithm> // std::min
14 #include <stdexcept> // std::invalid_argument
15 #include <iostream>
16 
17 #include "Util/SelectionRule.h"
18 #include "Util/CompInfo.h"
19 #include "LinAlg/SearchSpace.h"
20 #include "LinAlg/RitzPairs.h"
21 
22 namespace Spectra {
23 
33 template <typename Derived, typename OpType>
35 {
36 protected:
37  using Index = Eigen::Index;
38  using Scalar = typename OpType::Scalar;
39  using Matrix = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
40  using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
41 
42  const OpType& m_matrix_operator; // object to conduct matrix operation,
43  // e.g. matrix-vector product
44  Index niter_ = 0;
45  const Index m_number_eigenvalues; // number of eigenvalues requested
46  Index m_max_search_space_size;
47  Index m_initial_search_space_size;
48  Index m_correction_size; // how many correction vectors are added in each iteration
49  RitzPairs<Scalar> m_ritz_pairs; // Ritz eigen pair structure
50  SearchSpace<Scalar> m_search_space; // search space
51 
52 private:
53  CompInfo m_info = CompInfo::NotComputed; // status of the computation
54 
55  void check_argument() const
56  {
57  if (m_number_eigenvalues < 1 || m_number_eigenvalues > m_matrix_operator.cols() - 1)
58  throw std::invalid_argument("nev must satisfy 1 <= nev <= n - 1, n is the size of matrix");
59  }
60 
61 public:
62  JDSymEigsBase(OpType& op, Index nev) :
63  m_matrix_operator(op),
64  m_number_eigenvalues(nev),
65  m_max_search_space_size(10 * m_number_eigenvalues),
66  m_initial_search_space_size(2 * m_number_eigenvalues),
67  m_correction_size(m_number_eigenvalues)
68  {
69  check_argument();
70  //TODO better input validation and checks
71  if (op.cols() < m_max_search_space_size)
72  {
73  m_max_search_space_size = op.cols();
74  }
75 
76  if (op.cols() < m_initial_search_space_size + m_correction_size)
77  {
78  m_initial_search_space_size = op.cols() / 3;
79  m_correction_size = op.cols() / 3;
80  }
81  }
82 
86  void set_max_search_space_size(Index max_search_space_size)
87  {
88  m_max_search_space_size = max_search_space_size;
89  }
93  void set_correction_size(Index correction_size)
94  {
95  m_correction_size = correction_size;
96  }
97 
101  void set_initial_search_space_size(Index initial_search_space_size)
102  {
103  m_initial_search_space_size = initial_search_space_size;
104  }
105 
109  virtual ~JDSymEigsBase() {}
110 
115  CompInfo info() const { return m_info; }
116 
120  Index num_iterations() const { return niter_; }
121 
122  Vector eigenvalues() const { return m_ritz_pairs.ritz_values().head(m_number_eigenvalues); }
123 
124  Matrix eigenvectors() const { return m_ritz_pairs.ritz_vectors().leftCols(m_number_eigenvalues); }
125 
126  Index compute(SortRule selection = SortRule::LargestMagn, Index maxit = 100,
127  Scalar tol = 100 * Eigen::NumTraits<Scalar>::dummy_precision())
128  {
129  Derived& derived = static_cast<Derived&>(*this);
130  Matrix intial_space = derived.setup_initial_search_space(selection);
131  return compute_with_guess(intial_space, selection, maxit, tol);
132  }
133 
134  Index compute_with_guess(const Eigen::Ref<const Matrix>& initial_space,
135  SortRule selection = SortRule::LargestMagn,
136  Index maxit = 100,
137  Scalar tol = 100 * Eigen::NumTraits<Scalar>::dummy_precision())
138 
139  {
140  m_search_space.initialize_search_space(initial_space);
141  niter_ = 0;
142  for (niter_ = 0; niter_ < maxit; niter_++)
143  {
144  bool do_restart = (m_search_space.size() > m_max_search_space_size);
145 
146  if (do_restart)
147  {
148  m_search_space.restart(m_ritz_pairs, m_initial_search_space_size);
149  }
150 
151  m_search_space.update_operator_basis_product(m_matrix_operator);
152 
153  Eigen::ComputationInfo small_problem_info = m_ritz_pairs.compute_eigen_pairs(m_search_space);
154  if (small_problem_info != Eigen::ComputationInfo::Success)
155  {
156  m_info = CompInfo::NumericalIssue;
157  break;
158  }
159  m_ritz_pairs.sort(selection);
160 
161  bool converged = m_ritz_pairs.check_convergence(tol, m_number_eigenvalues);
162  if (converged)
163  {
164  m_info = CompInfo::Successful;
165  break;
166  }
167  else if (niter_ == maxit - 1)
168  {
169  m_info = CompInfo::NotConverging;
170  break;
171  }
172  Derived& derived = static_cast<Derived&>(*this);
173  Matrix corr_vect = derived.calculate_correction_vector();
174 
175  m_search_space.extend_basis(corr_vect);
176  }
177  return (m_ritz_pairs.converged_eigenvalues()).template cast<Index>().head(m_number_eigenvalues).sum();
178  }
179 };
180 
181 } // namespace Spectra
182 
183 #endif // SPECTRA_JD_SYM_EIGS_BASE_H
CompInfo info() const
void set_correction_size(Index correction_size)
Definition: JDSymEigsBase.h:93
Index num_iterations() const
void set_initial_search_space_size(Index initial_search_space_size)
void set_max_search_space_size(Index max_search_space_size)
Definition: JDSymEigsBase.h:86
CompInfo
Definition: CompInfo.h:18
@ Successful
Computation was successful.