Spectra  1.0.1 Header-only C++ Library for Large Scale Eigenvalue Problems
DavidsonSymEigsSolver.h
1 // Copyright (C) 2020 Netherlands eScience Center <f.zapata@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_DAVIDSON_SYM_EIGS_SOLVER_H
8 #define SPECTRA_DAVIDSON_SYM_EIGS_SOLVER_H
9
10 #include <Eigen/Core>
11
12 #include "JDSymEigsBase.h"
13 #include "Util/SelectionRule.h"
14
15 namespace Spectra {
16
30 template <typename OpType>
31 class DavidsonSymEigsSolver : public JDSymEigsBase<DavidsonSymEigsSolver<OpType>, OpType>
32 {
33 private:
34  using Index = Eigen::Index;
35  using Scalar = typename OpType::Scalar;
36  using Matrix = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
37  using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
38
39  Vector m_diagonal;
40
41 public:
42  DavidsonSymEigsSolver(OpType& op, Index nev, Index nvec_init, Index nvec_max) :
43  JDSymEigsBase<DavidsonSymEigsSolver<OpType>, OpType>(op, nev, nvec_init, nvec_max)
44  {
45  m_diagonal.resize(this->m_matrix_operator.rows());
46  for (Index i = 0; i < op.rows(); i++)
47  {
48  m_diagonal(i) = op(i, i);
49  }
50  }
51
52  DavidsonSymEigsSolver(OpType& op, Index nev) :
53  DavidsonSymEigsSolver(op, nev, 2 * nev, 10 * nev) {}
54
60  Matrix setup_initial_search_space(SortRule selection) const
61  {
62  std::vector<Eigen::Index> indices_sorted = argsort(selection, m_diagonal);
63
64  Matrix initial_basis = Matrix::Zero(this->m_matrix_operator.rows(), this->m_initial_search_space_size);
65
66  for (Index k = 0; k < this->m_initial_search_space_size; k++)
67  {
68  Index row = indices_sorted[k];
69  initial_basis(row, k) = 1.0;
70  }
71  return initial_basis;
72  }
73
78  {
79  const Matrix& residues = this->m_ritz_pairs.residues();
80  const Vector& eigvals = this->m_ritz_pairs.ritz_values();
81  Matrix correction = Matrix::Zero(this->m_matrix_operator.rows(), this->m_correction_size);
82  for (Index k = 0; k < this->m_correction_size; k++)
83  {
84  Vector tmp = eigvals(k) - m_diagonal.array();
85  correction.col(k) = residues.col(k).array() / tmp.array();
86  }
87  return correction;
88  }
89 };
90
91 } // namespace Spectra
92
93 #endif // SPECTRA_DAVIDSON_SYM_EIGS_SOLVER_H
Matrix setup_initial_search_space(SortRule selection) const