ergo
get_eigenvectors.h
Go to the documentation of this file.
1/* Ergo, version 3.8, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2019 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29#ifndef EIGENVECTORS_HEADER
30#define EIGENVECTORS_HEADER
31
32/************************************************/
33/* EIGENVECTORS COMPUTATIONS */
34/************************************************/
35
36
47#include "matrix_utilities.h"
49#include "SizesAndBlocks.h"
50#include "output.h"
51
52#include <iostream>
53#include <string.h>
54
56
57
58namespace eigvec
59{
61template<typename Treal, typename MatrixType, typename VectorType>
63{
64 VectorType y, Ay;
65
66 y = eigVec;
67 Treal ONE = (Treal)1.0;
68 y *= ONE / y.eucl(); // y = y/norm(y)
69 Ay = ONE * A * y; // Ay = A*y
70 Treal lambda = transpose(y) * Ay; // lambda = y'*Ay
71 return lambda;
72}
73
74
77template<typename Treal, typename MatrixType, typename VectorType>
79 std::vector<Treal>& eigVal,
80 std::vector<VectorType>& eigVec,
81 int number_of_eigenvalues,
82 const Treal TOL,
83 std::vector<int>& num_iter,
84 int maxit = 200,
85 bool do_deflation = false)
86{
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);
91
92 if (eigVec[0].is_empty())
93 {
94 throw std::runtime_error("Error in eigvec::lanczos_method : eigVec[0].is_empty()");
95 }
96
97 const Treal ONE = 1.0;
98
99 if (!do_deflation)
100 {
101 VectorType y;
102 y = eigVec[0];
103 y *= (ONE / y.eucl()); // normalization
104
106 (A, y, number_of_eigenvalues, maxit);
107
108 try
109 {
110 lan.setAbsTol(TOL);
111 lan.setRelTol(TOL);
112 lan.run();
113 }
114 catch (...)
115 {
116 num_iter[0] = maxit; // lanczos did not converge in maxIter iterations
117 throw;
118 }
119
120 Treal acc = 0;
121 for (int i = 1; i <= number_of_eigenvalues; i++) {
122 lan.get_ith_eigenpair(i, eigVal[i-1], eigVec[i-1], acc);
123 }
124
125 num_iter[0] = lan.get_num_iter();
126 }
127 else // do_deflation
128 {
129 // use the vector stored in eigVec[0]
130 if (number_of_eigenvalues > 1)
131 {
132 VectorType y, v1;
133 v1 = eigVec[0];
134
135 // get initial guess
136 if (eigVec[1].is_empty())
137 {
138 throw std::runtime_error("Error in eigvec::lanczos_method : eigVec[1].is_empty()");
139 }
140 y = eigVec[1];
141 y *= (ONE / y.eucl()); // normalization
142
143 try
144 {
145 int num_eig = 1; // just one eigenpair should be computed
146 // find bounds of the spectrum
147 Treal eigmin, eigmax;
148 A.gershgorin(eigmin, eigmax);
149 Treal sigma = eigVal[0] - eigmin; // out eigenvalue to the uninteresting end of the spectrum
151 (A, y, num_eig, maxit, 100, &v1, sigma);
152 lan.setAbsTol(TOL);
153 lan.setRelTol(TOL);
154 lan.run();
155 Treal acc = 0;
156 lan.get_ith_eigenpair(1, eigVal[1], eigVec[1], acc);
157
158 VectorType resVec(eigVec[1]); // residual
159 resVec *= eigVal[1];
160 resVec += -ONE * A * eigVec[1];
161
162 num_iter[1] = lan.get_num_iter();
163 }
164 catch (...)
165 {
166 num_iter[1] = maxit; // lanczos did not converge in maxIter iterations
167 throw;
168 }
169 }
170 else
171 {
172 throw std::runtime_error("Error in eigvec::lanczos_method : number_of_eigenvalues <= 1");
173 }
174 }
175}
176
177
180template<typename Treal, typename MatrixType, typename VectorType>
182 Treal& eigVal,
183 VectorType& eigVec,
184 const Treal TOL,
185 int& num_iter,
186 int maxit = 200)
187{
188 VectorType y;
189 VectorType Ay;
190 VectorType residual;
191 VectorType temp;
192 Treal lambda;
193 const Treal ONE = 1.0;
194 const Treal MONE = -1.0;
195
196 y = eigVec; // initial guess
197 y *= (ONE / y.eucl()); // normalization
198
199 // init
200 Ay = y;
201 residual = y;
202 temp = y;
203
204 int it = 1;
205 Ay = ONE * A * y; // Ay = A*y
206
207 while (it == 1 || (residual.eucl() / template_blas_fabs(lambda) > TOL && it <= maxit))
208 {
209 y = Ay;
210 y *= ONE / Ay.eucl(); // y = Ay/norm(Ay)
211 Ay = ONE * A * y; // Ay = A*y
212 lambda = transpose(y) * Ay; // lambda = y'*Ay
213
214 // r = A*y - lambda*y
215 residual = Ay;
216 residual += (MONE * lambda) * y;
217 //printf("residual.eucl() = %e\n", residual.eucl());
218
219 ++it;
220 }
221
222 printf("Power method required %d iterations.\n", it - 1);
223
224 eigVal = lambda;
225 eigVec = y;
226 num_iter = it - 1;
227}
228
229
231template<typename Treal, typename MatrixType, typename VectorType>
233 Treal tol,
234 std::vector<Treal>& eigVal,
235 std::vector<VectorType>& eigVec,
236 int number_of_eigenvalues_to_compute,
237 std::string method,
238 std::vector<int>& num_iter,
239 int maxit = 200,
240 bool do_deflation = false
241 )
242{
243 assert(number_of_eigenvalues_to_compute >= 1);
244 assert(eigVal.size() >= 1); // note: number_of_eigenvalues may not be equal to eigVal.size()
245 assert(eigVec.size() == eigVal.size());
246 assert(eigVec.size() == num_iter.size());
247
248 if (method == "power")
249 {
250 if (eigVal.size() > 1)
251 {
252 throw std::runtime_error("Error in eigvec::computeEigenvectors: computation of more "
253 "than 1 eigenpair is not implemented for the power method.");
254 }
255 if (do_deflation)
256 {
257 throw std::runtime_error("Error in eigvec::computeEigenvectors: deflation is not implemented for the power method.");
258 }
259 power_method(A, eigVal[0], eigVec[0], tol, num_iter[0], maxit);
260 }
261 else if (method == "lanczos")
262 {
263 lanczos_method(A, eigVal, eigVec, number_of_eigenvalues_to_compute, tol, num_iter, maxit, do_deflation);
264 }
265 else
266 {
267 throw std::runtime_error("Error in eigvec::computeEigenvectors: unknown method.");
268 }
269 return 0;
270}
271} // namespace
272
273#endif // EIGENVECTORS_HEADER
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 ...
Definition: VectorGeneral.h:48
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...
#define A
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)