ergo
get_eigenvectors.h
Go to the documentation of this file.
1 /* Ergo, version 3.7, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2018 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(number_of_eigenvalues >= 1);
90 
91  if (eigVec[0].is_empty())
92  {
93  throw std::runtime_error("Error in eigvec::lanczos_method : eigVec[0].is_empty()");
94  }
95 
96  const Treal ONE = 1.0;
97 
98  if (!do_deflation)
99  {
100  try
101  {
102  VectorType y;
103  y = eigVec[0];
104  y *= (ONE / y.eucl()); // normalization
105 
106 
108  (A, y, number_of_eigenvalues, maxit);
109  lan.setAbsTol(TOL);
110  lan.setRelTol(TOL);
111  lan.run();
112  Treal acc = 0;
113  lan.get_ith_eigenpair(1, eigVal[0], eigVec[0], acc);
114 
115  VectorType resVec(eigVec[0]); // residual
116  resVec *= eigVal[0];
117  resVec += -ONE * A * eigVec[0];
118 
119  /* if(number_of_eigenvalues == 2) */
120  /* { */
121  /* lan.get_ith_eigenpair(2, eigVal[1], eigVec[1], acc); */
122  /* VectorType resVec2(eigVec[1]); // residual */
123  /* resVec2 *= eigVal[1]; */
124  /* resVec2 += -ONE * A.MATRIX * eigVec[1]; */
125  /* } */
126  num_iter[0] = lan.get_num_iter();
127  }
128  catch (std::exception& e)
129  {
130  num_iter[0] = maxit; // lanczos did not converge in maxIter iterations
131  }
132  }
133  else // do_deflation
134  {
135  // use the vector stored in eigVec[0]
136  if (number_of_eigenvalues > 1)
137  {
138  VectorType y, v1;
139  v1 = eigVec[0];
140 
141  // get initial guess
142  if (eigVec[1].is_empty())
143  {
144  throw std::runtime_error("Error in eigvec::lanczos_method : eigVec[1].is_empty()");
145  }
146  y = eigVec[1];
147  y *= (ONE / y.eucl()); // normalization
148 
149  try
150  {
151  int num_eig = 1; // just one eigenpair should be computed
152  // find bounds of the spectrum
153  Treal eigmin, eigmax;
154  A.gershgorin(eigmin, eigmax);
155  Treal sigma = eigVal[0] - eigmin; // out eigenvalue to the uninteresting end of the spectrum
157  (A, y, num_eig, maxit, 100, &v1, sigma);
158  lan.setAbsTol(TOL);
159  lan.setRelTol(TOL);
160  lan.run();
161  Treal acc = 0;
162  lan.get_ith_eigenpair(1, eigVal[1], eigVec[1], acc);
163 
164  VectorType resVec(eigVec[1]); // residual
165  resVec *= eigVal[1];
166  resVec += -ONE * A * eigVec[1];
167 
168  num_iter[1] = lan.get_num_iter();
169  }
170  catch (std::exception& e)
171  {
172  num_iter[1] = maxit; // lanczos did not converge in maxIter iterations
173  }
174  }
175  else
176  {
177  throw std::runtime_error("Error in eigvec::lanczos_method : number_of_eigenvalues <= 1");
178  }
179  }
180 }
181 
182 
185 template<typename Treal, typename MatrixType, typename VectorType>
187  Treal& eigVal,
188  VectorType& eigVec,
189  const Treal TOL,
190  int& num_iter,
191  int maxit = 200)
192 {
193  VectorType y;
194  VectorType Ay;
195  VectorType residual;
196  VectorType temp;
197  Treal lambda;
198  const Treal ONE = 1.0;
199  const Treal MONE = -1.0;
200 
201  y = eigVec; // initial guess
202  y *= (ONE / y.eucl()); // normalization
203 
204  // init
205  Ay = y;
206  residual = y;
207  temp = y;
208 
209  int it = 1;
210  Ay = ONE * A * y; // Ay = A*y
211 
212  while (it == 1 || (residual.eucl() / template_blas_fabs(lambda) > TOL && it <= maxit))
213  {
214  y = Ay;
215  y *= ONE / Ay.eucl(); // y = Ay/norm(Ay)
216  Ay = ONE * A * y; // Ay = A*y
217  lambda = transpose(y) * Ay; // lambda = y'*Ay
218 
219  // r = A*y - lambda*y
220  residual = Ay;
221  residual += (MONE * lambda) * y;
222  //printf("residual.eucl() = %e\n", residual.eucl());
223 
224  ++it;
225  }
226 
227  printf("Power method required %d iterations.\n", it - 1);
228 
229  eigVal = lambda;
230  eigVec = y;
231  num_iter = it - 1;
232 }
233 
234 
236 template<typename Treal, typename MatrixType, typename VectorType>
238  Treal tol,
239  std::vector<Treal>& eigVal,
240  std::vector<VectorType>& eigVec,
241  int number_of_eigenvalues_to_compute,
242  std::string method,
243  std::vector<int>& num_iter,
244  int maxit = 200,
245  bool do_deflation = false
246  )
247 {
248  assert(number_of_eigenvalues_to_compute >= 1);
249  assert(eigVal.size() >= 1); // note: number_of_eigenvalues may not be equal to eigVal.size()
250  assert(eigVec.size() == eigVal.size());
251  assert(eigVec.size() == num_iter.size());
252 
253  if (method == "power")
254  {
255  if (eigVal.size() > 1)
256  {
257  throw std::runtime_error("Error in eigvec::computeEigenvectors: computation of more "
258  "than 1 eigenpair is not implemented for the power method.");
259  }
260  if (do_deflation)
261  {
262  throw std::runtime_error("Error in eigvec::computeEigenvectors: deflation is not implemented for the power method.");
263  }
264  power_method(A, eigVal[0], eigVec[0], tol, num_iter[0], maxit);
265  }
266  else if (method == "lanczos")
267  {
268  lanczos_method(A, eigVal, eigVec, number_of_eigenvalues_to_compute, tol, num_iter, maxit, do_deflation);
269  }
270  else
271  {
272  throw std::runtime_error("Error in eigvec::computeEigenvectors: unknown method.");
273  }
274  return 0;
275 }
276 } // namespace
277 
278 #endif // EIGENVECTORS_HEADER
#define A
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&#39;Ay)/(y&#39;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