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 
58 namespace eigvec
59 {
61 template<typename Treal, typename MatrixType, typename VectorType>
62 Treal compute_rayleigh_quotient(const MatrixType& A, const VectorType& eigVec)
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 
77 template<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 
180 template<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 
231 template<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
eigvec::compute_rayleigh_quotient
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
integral_matrix_wrappers.h
Wrapper routines for different parts of the integral code, including conversion of matrices from/to t...
eigvec::computeEigenvectors
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
mat::arn::LanczosSeveralLargestEig::setRelTol
void setRelTol(Treal const newTol)
Definition: LanczosSeveralLargestEig.h:100
mat::VectorGeneral::eucl
Treal eucl() const
Definition: VectorGeneral.h:178
template_blas_fabs
Treal template_blas_fabs(Treal x)
eigvec::lanczos_method
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
mat::arn::LanczosSeveralLargestEig
Definition: LanczosSeveralLargestEig.h:51
MatrixType
symmMatrix MatrixType
Definition: recexp_diag_test.cc:66
matrix_utilities.h
Utilities related to the hierarchical matrix library (HML), including functions for setting up permut...
mat::VectorGeneral
Definition: MatrixBase.h:55
eigvec
Definition: get_eigenvectors.h:59
mat::transpose
Xtrans< TX > transpose(TX const &A)
Transposition.
Definition: matrix_proxy.h:131
A
#define A
LanczosSeveralLargestEig.h
mat::arn::LanczosSeveralLargestEig::setAbsTol
void setAbsTol(Treal const newTol)
Definition: LanczosSeveralLargestEig.h:101
mat::arn::LanczosSeveralLargestEig::get_num_iter
int get_num_iter() const
Definition: LanczosSeveralLargestEig.h:159
mat::arn::LanczosSeveralLargestEig::run
virtual void run()
Definition: LanczosSeveralLargestEig.h:108
SizesAndBlocks.h
Class used to keep track of the block sizes used at different levels in the hierarchical matrix data ...
eigvec::power_method
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
mat::arn::LanczosSeveralLargestEig::get_ith_eigenpair
virtual void get_ith_eigenpair(int i, Treal &eigVal, Tvector &eigVec, Treal &acc)
Definition: LanczosSeveralLargestEig.h:149
output.h
Functionality for writing output messages to a text file.