29 #ifndef EIGENVECTORS_HEADER 30 #define EIGENVECTORS_HEADER 61 template<
typename Treal,
typename MatrixType,
typename VectorType>
67 Treal ONE = (Treal)1.0;
77 template<
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(number_of_eigenvalues >= 1);
91 if (eigVec[0].is_empty())
93 throw std::runtime_error(
"Error in eigvec::lanczos_method : eigVec[0].is_empty()");
96 const Treal ONE = 1.0;
104 y *= (ONE / y.
eucl());
108 (
A, y, number_of_eigenvalues, maxit);
117 resVec += -ONE *
A * eigVec[0];
128 catch (std::exception& e)
136 if (number_of_eigenvalues > 1)
142 if (eigVec[1].is_empty())
144 throw std::runtime_error(
"Error in eigvec::lanczos_method : eigVec[1].is_empty()");
147 y *= (ONE / y.
eucl());
153 Treal eigmin, eigmax;
154 A.gershgorin(eigmin, eigmax);
155 Treal sigma = eigVal[0] - eigmin;
157 (
A, y, num_eig, maxit, 100, &v1, sigma);
166 resVec += -ONE *
A * eigVec[1];
170 catch (std::exception& e)
177 throw std::runtime_error(
"Error in eigvec::lanczos_method : number_of_eigenvalues <= 1");
185 template<
typename Treal,
typename MatrixType,
typename VectorType>
198 const Treal ONE = 1.0;
199 const Treal MONE = -1.0;
202 y *= (ONE / y.
eucl());
215 y *= ONE / Ay.
eucl();
221 residual += (MONE * lambda) * y;
227 printf(
"Power method required %d iterations.\n", it - 1);
236 template<
typename Treal,
typename MatrixType,
typename VectorType>
239 std::vector<Treal>& eigVal,
240 std::vector<VectorType>& eigVec,
241 int number_of_eigenvalues_to_compute,
243 std::vector<int>& num_iter,
245 bool do_deflation =
false 248 assert(number_of_eigenvalues_to_compute >= 1);
249 assert(eigVal.size() >= 1);
250 assert(eigVec.size() == eigVal.size());
251 assert(eigVec.size() == num_iter.size());
253 if (method ==
"power")
255 if (eigVal.size() > 1)
257 throw std::runtime_error(
"Error in eigvec::computeEigenvectors: computation of more " 258 "than 1 eigenpair is not implemented for the power method.");
262 throw std::runtime_error(
"Error in eigvec::computeEigenvectors: deflation is not implemented for the power method.");
264 power_method(
A, eigVal[0], eigVec[0], tol, num_iter[0], maxit);
266 else if (method ==
"lanczos")
268 lanczos_method(
A, eigVal, eigVec, number_of_eigenvalues_to_compute, tol, num_iter, maxit, do_deflation);
272 throw std::runtime_error(
"Error in eigvec::computeEigenvectors: unknown method.");
278 #endif // EIGENVECTORS_HEADER
Definition: LanczosSeveralLargestEig.h:50
void setRelTol(Treal const newTol)
Definition: LanczosSeveralLargestEig.h:100
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:186
Functionality for writing output messages to a text file.
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:237
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...
Class used to keep track of the block sizes used at different levels in the hierarchical matrix data ...
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 template_blas_fabs(Treal x)
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
symmMatrix MatrixType
Definition: recexp_diag_test.cc:66
Class for computing several largest (note: not by magnitude) eigenvalues of a symmetric matrix with t...
virtual void get_ith_eigenpair(int i, Treal &eigVal, Tvector &eigVec, Treal &acc)
Definition: LanczosSeveralLargestEig.h:124
Definition: get_eigenvectors.h:58
int get_num_iter() const
Definition: LanczosSeveralLargestEig.h:134
Treal eucl() const
Definition: VectorGeneral.h:172
void setAbsTol(Treal const newTol)
Definition: LanczosSeveralLargestEig.h:101
virtual void run()
Definition: LanczosSeveralLargestEig.h:108
Xtrans< TX > transpose(TX const &A)
Transposition.
Definition: matrix_proxy.h:131