29#ifndef EIGENVECTORS_HEADER
30#define EIGENVECTORS_HEADER
61template<
typename Treal,
typename MatrixType,
typename VectorType>
67 Treal ONE = (Treal)1.0;
70 Treal lambda = transpose(y) * Ay;
77template<
typename Treal,
typename MatrixType,
typename VectorType>
79 std::vector<Treal>& eigVal,
80 std::vector<VectorType>& eigVec,
81 int number_of_eigenvalues,
83 std::vector<int>& num_iter,
85 bool do_deflation =
false)
87 assert(eigVal.size() == eigVec.size());
88 assert(eigVal.size() == num_iter.size());
89 assert((
int) eigVal.size() >= number_of_eigenvalues);
90 assert(number_of_eigenvalues >= 1);
92 if (eigVec[0].is_empty())
94 throw std::runtime_error(
"Error in eigvec::lanczos_method : eigVec[0].is_empty()");
97 const Treal ONE = 1.0;
103 y *= (ONE / y.
eucl());
106 (
A, y, number_of_eigenvalues, maxit);
121 for (
int i = 1; i <= number_of_eigenvalues; i++) {
130 if (number_of_eigenvalues > 1)
136 if (eigVec[1].is_empty())
138 throw std::runtime_error(
"Error in eigvec::lanczos_method : eigVec[1].is_empty()");
141 y *= (ONE / y.
eucl());
147 Treal eigmin, eigmax;
148 A.gershgorin(eigmin, eigmax);
149 Treal sigma = eigVal[0] - eigmin;
151 (
A, y, num_eig, maxit, 100, &v1, sigma);
160 resVec += -ONE *
A * eigVec[1];
172 throw std::runtime_error(
"Error in eigvec::lanczos_method : number_of_eigenvalues <= 1");
180template<
typename Treal,
typename MatrixType,
typename VectorType>
193 const Treal ONE = 1.0;
194 const Treal MONE = -1.0;
197 y *= (ONE / y.
eucl());
210 y *= ONE / Ay.
eucl();
212 lambda = transpose(y) * Ay;
216 residual += (MONE * lambda) * y;
222 printf(
"Power method required %d iterations.\n", it - 1);
231template<
typename Treal,
typename MatrixType,
typename VectorType>
234 std::vector<Treal>& eigVal,
235 std::vector<VectorType>& eigVec,
236 int number_of_eigenvalues_to_compute,
238 std::vector<int>& num_iter,
240 bool do_deflation =
false
243 assert(number_of_eigenvalues_to_compute >= 1);
244 assert(eigVal.size() >= 1);
245 assert(eigVec.size() == eigVal.size());
246 assert(eigVec.size() == num_iter.size());
248 if (method ==
"power")
250 if (eigVal.size() > 1)
252 throw std::runtime_error(
"Error in eigvec::computeEigenvectors: computation of more "
253 "than 1 eigenpair is not implemented for the power method.");
257 throw std::runtime_error(
"Error in eigvec::computeEigenvectors: deflation is not implemented for the power method.");
259 power_method(
A, eigVal[0], eigVec[0], tol, num_iter[0], maxit);
261 else if (method ==
"lanczos")
263 lanczos_method(
A, eigVal, eigVec, number_of_eigenvalues_to_compute, tol, num_iter, maxit, do_deflation);
267 throw std::runtime_error(
"Error in eigvec::computeEigenvectors: unknown method.");
generalVector VectorType
Definition GetDensFromFock.cc:62
Class for computing several largest (note: not by magnitude) eigenvalues of a symmetric matrix with t...
Class used to keep track of the block sizes used at different levels in the hierarchical matrix data ...
Treal eucl() const
Definition VectorGeneral.h:178
Definition LanczosSeveralLargestEig.h:51
int get_num_iter() const
Definition LanczosSeveralLargestEig.h:159
virtual void run()
Definition LanczosSeveralLargestEig.h:108
void setAbsTol(Treal const newTol)
Definition LanczosSeveralLargestEig.h:101
virtual void get_ith_eigenpair(int i, Treal &eigVal, Tvector &eigVec, Treal &acc)
Definition LanczosSeveralLargestEig.h:149
void setRelTol(Treal const newTol)
Definition LanczosSeveralLargestEig.h:100
Wrapper routines for different parts of the integral code, including conversion of matrices from/to t...
Utilities related to the hierarchical matrix library (HML), including functions for setting up permut...
Definition get_eigenvectors.h:59
void lanczos_method(const MatrixType &A, std::vector< Treal > &eigVal, std::vector< VectorType > &eigVec, int number_of_eigenvalues, const Treal TOL, std::vector< int > &num_iter, int maxit=200, bool do_deflation=false)
Use Lanzcos method for computing eigenvectors.
Definition get_eigenvectors.h:78
Treal compute_rayleigh_quotient(const MatrixType &A, const VectorType &eigVec)
Get Rayleigh quotient: A = (y'Ay)/(y'y), y = eigVecPtr.
Definition get_eigenvectors.h:62
int computeEigenvectors(const MatrixType &A, Treal tol, std::vector< Treal > &eigVal, std::vector< VectorType > &eigVec, int number_of_eigenvalues_to_compute, std::string method, std::vector< int > &num_iter, int maxit=200, bool do_deflation=false)
Function for choosing method for computing eigenvectors.
Definition get_eigenvectors.h:232
void power_method(const MatrixType &A, Treal &eigVal, VectorType &eigVec, const Treal TOL, int &num_iter, int maxit=200)
Use power method for computing eigenvectors.
Definition get_eigenvectors.h:181
Functionality for writing output messages to a text file.
symmMatrix MatrixType
Definition recexp_diag_test.cc:66
Treal template_blas_fabs(Treal x)