Spectra  1.0.0
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) :
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 
57  Matrix setup_initial_search_space(SortRule selection) const
58  {
59  std::vector<Eigen::Index> indices_sorted = argsort(selection, m_diagonal);
60 
61  Matrix initial_basis = Matrix::Zero(this->m_matrix_operator.rows(), this->m_initial_search_space_size);
62 
63  for (Index k = 0; k < this->m_initial_search_space_size; k++)
64  {
65  Index row = indices_sorted[k];
66  initial_basis(row, k) = 1.0;
67  }
68  return initial_basis;
69  }
70 
75  {
76  const Matrix& residues = this->m_ritz_pairs.residues();
77  const Vector& eigvals = this->m_ritz_pairs.ritz_values();
78  Matrix correction = Matrix::Zero(this->m_matrix_operator.rows(), this->m_correction_size);
79  for (Index k = 0; k < this->m_correction_size; k++)
80  {
81  Vector tmp = eigvals(k) - m_diagonal.array();
82  correction.col(k) = residues.col(k).array() / tmp.array();
83  }
84  return correction;
85  }
86 };
87 
88 } // namespace Spectra
89 
90 #endif // SPECTRA_DAVIDSON_SYM_EIGS_SOLVER_H
Matrix setup_initial_search_space(SortRule selection) const