ergo
LanczosSeveralLargestEig.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 
40 #ifndef MAT_LANCZOSSEVERALLARGESTMAGNITUDEEIG
41 #define MAT_LANCZOSSEVERALLARGESTMAGNITUDEEIG
42 
43 #include <limits>
44 #include <vector>
45 
46 namespace mat { /* Matrix namespace */
47  namespace arn { /* Arnoldi type methods namespace */
48 
49  template<typename Treal, typename Tmatrix, typename Tvector>
51  {
52  public:
53  // AA - matrix
54  // startVec - starting guess vector
55  // num_eigs - number of eigenvalues to compute
56  // maxIter(100) - number of iterations
57  // cap(100) - estimated number of vectors in the Krylov subspace, will be increased if needed automatically
58  // deflVec_(NULL) - (deflation) vector corresponding to an uninteresting eigenvalue
59  // sigma_(0) - (deflation) shift of an uninteresting eigenvalue (to put it in the uninteresting part of the spectrum)
60  LanczosSeveralLargestEig(Tmatrix const & AA, Tvector const & startVec, int num_eigs,
61  int maxit = 100, int cap = 100, Tvector * deflVec_ = NULL, Treal sigma_ = 0)
62  : A(AA),
63  v(new Tvector[cap]),
64  eigVectorTri(0),
65  capacity(cap),
66  j(0),
67  maxIter(maxit),
68  eValTmp(0),
69  accTmp(0),
70  number_of_eigenv(num_eigs),
71  alpha(0),
72  beta(0),
73  use_selective_orth(false),
74  use_full_orth(true),
75  counter_all(0),
76  counter_orth(0),
77  deflVec(deflVec_)
78  {
79  assert(cap > 1);
80  Treal const ONE = 1.0;
81  v[0] = startVec;
82  if(v[0].eucl() < template_blas_sqrt(getRelPrecision<Treal>())) {
83  v[0].rand();
84  }
85  v[0] *= (ONE / v[0].eucl());
86  r = v[0];
87 
88  if(number_of_eigenv == 1)
90 
91  absTol = 1e-12;
92  relTol = 1e-12;
93  sigma = sigma_;
94 
95  }
96 
97  // Absolute and relative tolerances
98  // Absolute accuracy is measured by the residual ||Ax-lambda*x||
99  // Realtive accuracy is measured by the relative residual ||Ax-lambda*x||/|lambda|
100  void setRelTol(Treal const newTol) { relTol = newTol; }
101  void setAbsTol(Treal const newTol) { absTol = newTol; }
102 
107 
108  virtual void run() {
109  do {
110  if(j > 1 && use_selective_orth)
111  selective_orth();
112  step();
113  update();
114  if (j > maxIter)
115  throw AcceptableMaxIter("Lanczos::run() did not converge within maxIter");
116  }
117  while (!converged());
118 
119  total_num_iter = j;
120 
121  // check orthogonality just in case
122  if(number_of_eigenv > 1)
123  {
124  for(int i = 0; i < total_num_iter-1; ++i)
125  for(int k = 0; k < total_num_iter-1; ++k)
126  {
127  if(i == k) continue;
128  v[i].readFromFile(); v[k].readFromFile();
129  Treal val = transpose(v[i]) * v[k]; // should be 0
130  if(val > template_blas_sqrt(mat::getMachineEpsilon<Treal>()))
131  throw std::runtime_error("Lanczos::run() : detected loss of orthogonality! Discard results.");
132  v[i].writeToFile(); v[k].writeToFile();
133  }
134  for(int k = 0; k < total_num_iter-1; ++k)
135  {
136  v[k].readFromFile();
137  Treal val = transpose(v[k]) * v[total_num_iter]; // should be 0
138  if(val > template_blas_sqrt(mat::getMachineEpsilon<Treal>()))
139  throw std::runtime_error("Lanczos::run() : detected loss of orthogonality! Discard results.");
140  v[k].writeToFile();
141  }
142  }
143  }
144 
145 
146 
147 
148  // i is a number of eigenvalue (1 is the largest, 2 is the second largest and so on)
149  virtual void get_ith_eigenpair(int i, Treal& eigVal, Tvector& eigVec, Treal & acc)
150  {
151  assert(i > 0);
152  assert(i <= size_accTmp);
153  eigVal = eValTmp[size_accTmp - i]; // array
154  assert(eigVectorTri);
155  getEigVector(eigVec, &eigVectorTri[j * (size_accTmp - i)]);
156  acc = accTmp[size_accTmp - i];
157  }
158 
159  int get_num_iter() const{ return total_num_iter;}
160 
162 
164  printf("Orthogonalized %d of total possible %d, this is %lf %%\n", counter_orth, counter_all, (double)counter_orth/counter_all*100);
165 
166  delete[] eigVectorTri;
167  delete[] eValTmp;
168  delete[] accTmp;
169  delete[] v;
170  }
171 
172 
174  Tricopy = Tri;
175  }
176 
177 
178  protected:
179  Tmatrix const & A;
180  Tvector* v;
185  Tvector r;
187  Treal* eigVectorTri; // Eigenvectors of the tridiagonal matrix
188  int capacity;
189  int j;
190  int maxIter;
191  void increaseCapacity(int const newCapacity);
192  void getEigVector(Tvector& eigVec, Treal const * const eVecTri) const;
193  Treal absTol;
194  Treal relTol;
195  virtual void step();
196  virtual void computeEigenPairTri();
197  virtual void update() {
199  }
200  void selective_orth();
201  virtual bool converged() const;
202  virtual bool converged_ith(int i) const;
203  Treal* eValTmp; // current computed eigenvalues (less or equal to number_of_eigenv)
204  Treal* accTmp; // residuals
205  int number_of_eigenv; // eigenvalues are saved in the decreasing order, thus the largest one has index 1
206  int size_accTmp; // size of accTmp (number of computed eigenvalues of the matrix T)
207  private:
208  Treal alpha;
209  Treal beta;
210 
212 
215 
218 
219  // if deflation is used
220  Tvector * deflVec;
221  Treal sigma;
222  };
223 
224  template<typename Treal, typename Tmatrix, typename Tvector>
227  {
228  int j_curr = j-1;
229 
230  Treal coeff = 0, res;
231  Treal normT = 0; // spectral norm of T (since norm of A is not available)
232  // find largest by absolute value eigenvalue of T
233  for(int i = 0; i <= j_curr; ++i)
234  if(template_blas_fabs(eValTmp[i]) > normT) normT = template_blas_fabs(eValTmp[i]);
235 
236  Treal epsilon = mat::getMachineEpsilon<Treal>();
237  Tvector tmp;
238  tmp = v[j_curr+1];
239  tmp *= beta; // return non-normalized value
240 
241  for(int i = j_curr; i >= 0; --i)
242  {
243  counter_all++;
244  // get residual for this eigenpair
245  res = accTmp[i];
246  Treal tol = template_blas_sqrt(epsilon) * normT;
247  if(res <= tol) // b_{j} * |VT_i(j)| <= sqrt(eps) * norm(A), but we do not have norm(A)
248  {
249  counter_orth++;
250  Tvector eigVec;
251  getEigVector(eigVec, &eigVectorTri[j_curr * i]); // y = U*VT(:, i); % ith Ritz vector
252  coeff = transpose(eigVec) * tmp;
253  tmp += (-coeff) * (eigVec); // v = v - (y'*v)*y
254  }
255  }
256 
257 
258 
259  v[j_curr+1] = tmp;
260  beta = v[j_curr+1].eucl(); // update beta
261  Treal const ONE = 1.0;
262  v[j_curr+1] *= ONE / beta; // normalized
263  Tri.update_beta(beta);
264  }
265 
266 
267 
268 
269  template<typename Treal, typename Tmatrix, typename Tvector>
271  step()
272  {
273  if (j + 1 >= capacity)
274  increaseCapacity(capacity * 2);
275  Treal const ONE = 1.0;
276  A.matVecProd(r, v[j]); // r = A * v[j]
277  alpha = transpose(v[j]) * r; // alpha = v[j]'*A*v[j]
278 
279  /*
280  If one wants to use deflation with vector
281  x_1:=deflVec (usually it is an eigenvector
282  corresponding to an eigenvalue lambda_1 of A)
283  and thus compute eigenvalues of the matrix
284  An = A-sigma*x_1*x_1'
285  Note: if lambga_i are eigenvalues of A corresponding to x_i, then
286  An will have eigenvalues (lambda_1-sigma, lambda_2, ..., lambda_N)
287  and unchanged eigenvectors x_i.
288  */
289 
290  if(deflVec != NULL)
291  {
292  /*
293  r = (A*vj - sigma*(x_1'*vj)*x_1) - alpha*vj - beta*v{j-1}
294  where
295  alpha = vj'*An*vj = vj'*A*vj - sigma * (x_1'*vj)^2
296  */
297  Treal gamma = transpose(*deflVec) * v[j]; // dot product x' * v_j
298  alpha -= sigma*gamma*gamma;
299 
300  r += (-sigma*gamma) * (*deflVec);
301  }
302 
303  r += (-alpha) * v[j];
304 
305  if (j) {
306  r += (-beta) * v[j-1];
307  v[j-1].writeToFile();
308  }
309 
310 
311  /*
312  If we need many eigenpairs, Lanczos vectors loose orthogonality as soon as one of the eigenpairs converges. If we continue iterations, then may appear some spurious eigenvalues. These spurious eigenvalues will eventually converge to the existing ones and we will get multiple convergence to the same eigenvalue. (In principle, we can probabaly check if we already converged to some eigenvalue before and just ignore it.) We use the simplest fix to the orthogonality loss, the full re-orthogonalization. This makes Lanczos procedure essentially equivalent to the Arnoldi algorithm. The only difference is that we are still using tridiagonal matrix.
313  */
314  if(use_full_orth)
315  {
316  // full re-orthogonalization (modified Gram-Schmidt)
317  Treal gamma_i = 0;
318  for(int i = 0; i < j; ++i )
319  {
320  v[i].readFromFile();
321  gamma_i = transpose(r) * v[i]; // r'*v_i
322  r += (-gamma_i) * v[i]; // (r'*vi) * v_i
323  v[i].writeToFile();
324  }
325  gamma_i = transpose(r) * v[j]; // r'*v_i
326  r += (-gamma_i) * v[j]; // (r'*vi) * v_i
327  }
328 
329 
330  beta = r.eucl();
331  v[j+1] = r;
332  v[j+1] *= ONE / beta;
333  Tri.increase(alpha, beta);
334  j++;
335  }
336 
337 
338  /*
339  Compute eigenvectors of the tridiagonal matrix
340  */
341  template<typename Treal, typename Tmatrix, typename Tvector>
344  if( eigVectorTri != NULL ) delete[] eigVectorTri;
345  if( accTmp != NULL ) delete[] accTmp;
346  if( eValTmp != NULL ) delete[] eValTmp;
347 
348  int num_compute_eigenvalues;
349  if(use_selective_orth)
350  num_compute_eigenvalues = j; // we need all eigenvectors of T
351  else
352  num_compute_eigenvalues = number_of_eigenv; // it is enough just number_of_eigenv of T
353 
354  /* Get largest eigenvalues */
355  int const max_ind = j-1; // eigenvalue count starts with 0
356  int const min_ind = std::max(j - num_compute_eigenvalues, 0);
357 
358  Treal* eigVectors = new Treal[j * num_compute_eigenvalues]; // every vector of size j
359  Treal* eigVals = new Treal[num_compute_eigenvalues];
360  Treal* accMax = new Treal[num_compute_eigenvalues];
361  assert(eigVectors != NULL);
362  assert(eigVals != NULL);
363  assert(accMax != NULL);
364 
365  Tri.getEigsByIndex(eigVals, eigVectors, accMax,
366  min_ind, max_ind);
367 
368  eValTmp = eigVals;
369 
370 
371  eigVectorTri = eigVectors;
372  accTmp = accMax;
373  size_accTmp = num_compute_eigenvalues;
374 
375  // set unused pointers to NULL
376  eigVectors = NULL;
377  eigVals = NULL;
378  accMax = NULL;
379  }
380 
381 
382 
383 
384  /* FIXME: If several eigenvectors are needed it is more optimal to
385  * compute all of them at once since then the krylov subspace vectors
386  * only need to be read from memory once.
387  */
388  template<typename Treal, typename Tmatrix, typename Tvector>
390  getEigVector(Tvector& eigVec, Treal const * const eVecTri) const {
391  if (j <= 1) {
392  eigVec = v[0];
393  }
394  else {
395  v[0].readFromFile();
396  eigVec = v[0];
397  v[0].writeToFile();
398  }
399  eigVec *= eVecTri[0];
400  for (int ind = 1; ind <= j - 2; ++ind) {
401  v[ind].readFromFile();
402  eigVec += eVecTri[ind] * v[ind];
403  v[ind].writeToFile();
404  }
405  eigVec += eVecTri[j-1] * v[j-1];
406 
407  // normalized
408  Treal norm_eigVec = eigVec.eucl();
409  Treal const ONE = 1.0;
410  eigVec *= ONE / norm_eigVec;
411  }
412 
413 
414  // we want lowest eigenvalue to converge
415  template<typename Treal, typename Tmatrix, typename Tvector>
417  converged() const {
418 
419  if(j < number_of_eigenv) return false;
420  bool conv1 = true;
421  if(number_of_eigenv > 1)
422  conv1 = converged_ith(number_of_eigenv-1);
423  bool conv = converged_ith(number_of_eigenv); // if the last needed eigenvalue converged
424 
425  return conv && conv1;
426  }
427 
428  // check convergence of ith eigenpair
429  template<typename Treal, typename Tmatrix, typename Tvector>
431  converged_ith(int i) const {
432  assert(size_accTmp >= i);
433 
434  bool conv = true; //accTmp[size_accTmp - i] < absTol; /* Do not use absolute accuracy */
435  if (template_blas_fabs(eValTmp[size_accTmp - i]) > 0) {
436  conv = conv &&
437  accTmp[size_accTmp - i] / template_blas_fabs(eValTmp[size_accTmp - i]) < relTol; /* Relative acc.*/
438  }
439  return conv;
440  }
441 
442 
443  template<typename Treal, typename Tmatrix, typename Tvector>
445  increaseCapacity(int const newCapacity) {
446  assert(newCapacity > capacity);
447  assert(j > 0);
448  capacity = newCapacity;
449  Tvector* new_v = new Tvector[capacity];
450  assert(new_v != NULL);
451  /* FIXME: Fix so that file is copied when operator= is called in Vector
452  * class
453  */
454  for (int ind = 0; ind <= j - 2; ind++){
455  v[ind].readFromFile();
456  new_v[ind] = v[ind];
457  new_v[ind].writeToFile();
458  }
459  for (int ind = j - 1; ind <= j; ind++){
460  new_v[ind] = v[ind];
461  }
462  delete[] v;
463  v = new_v;
464  }
465 
466 
467  } /* end namespace arn */
468 
469 
470 } /* end namespace mat */
471 #endif
mat::AcceptableMaxIter
Definition: Failure.h:81
mat::arn::LanczosSeveralLargestEig::counter_all
int counter_all
Definition: LanczosSeveralLargestEig.h:216
mat::arn::LanczosSeveralLargestEig::LanczosSeveralLargestEig
LanczosSeveralLargestEig(Tmatrix const &AA, Tvector const &startVec, int num_eigs, int maxit=100, int cap=100, Tvector *deflVec_=NULL, Treal sigma_=0)
Definition: LanczosSeveralLargestEig.h:60
template_blas_sqrt
Treal template_blas_sqrt(Treal x)
mat::arn::LanczosSeveralLargestEig::j
int j
Definition: LanczosSeveralLargestEig.h:189
mat::arn::LanczosSeveralLargestEig::copyTridiag
void copyTridiag(MatrixTridiagSymmetric< Treal > &Tricopy)
Definition: LanczosSeveralLargestEig.h:173
mat::arn::LanczosSeveralLargestEig::deflVec
Tvector * deflVec
Definition: LanczosSeveralLargestEig.h:220
mat::arn::LanczosSeveralLargestEig::setRelTol
void setRelTol(Treal const newTol)
Definition: LanczosSeveralLargestEig.h:100
mat::arn::LanczosSeveralLargestEig::total_num_iter
int total_num_iter
Definition: LanczosSeveralLargestEig.h:211
mat::arn::LanczosSeveralLargestEig::A
Tmatrix const & A
Definition: LanczosSeveralLargestEig.h:179
mat::arn::LanczosSeveralLargestEig::set_use_selective_orth
void set_use_selective_orth()
Definition: LanczosSeveralLargestEig.h:103
mat::arn::LanczosSeveralLargestEig::use_full_orth
bool use_full_orth
Definition: LanczosSeveralLargestEig.h:214
mat::arn::LanczosSeveralLargestEig::computeEigenPairTri
virtual void computeEigenPairTri()
Definition: LanczosSeveralLargestEig.h:343
mat::arn::LanczosSeveralLargestEig::converged
virtual bool converged() const
Definition: LanczosSeveralLargestEig.h:417
mat::arn::LanczosSeveralLargestEig::absTol
Treal absTol
Definition: LanczosSeveralLargestEig.h:193
mat::arn::LanczosSeveralLargestEig::beta
Treal beta
Definition: LanczosSeveralLargestEig.h:209
template_blas_fabs
Treal template_blas_fabs(Treal x)
mat::arn::LanczosSeveralLargestEig::sigma
Treal sigma
Definition: LanczosSeveralLargestEig.h:221
mat::arn::LanczosSeveralLargestEig::number_of_eigenv
int number_of_eigenv
Definition: LanczosSeveralLargestEig.h:205
mat::arn::LanczosSeveralLargestEig::maxIter
int maxIter
Current step.
Definition: LanczosSeveralLargestEig.h:190
mat::arn::LanczosSeveralLargestEig::capacity
int capacity
Definition: LanczosSeveralLargestEig.h:188
mat::arn::MatrixTridiagSymmetric
Tridiagonal symmetric matrix class template.
Definition: MatrixTridiagSymmetric.h:47
mat::arn::LanczosSeveralLargestEig
Definition: LanczosSeveralLargestEig.h:51
mat::arn::LanczosSeveralLargestEig::Tri
MatrixTridiagSymmetric< Treal > Tri
Residual vector.
Definition: LanczosSeveralLargestEig.h:186
mat::arn::LanczosSeveralLargestEig::increaseCapacity
void increaseCapacity(int const newCapacity)
Definition: LanczosSeveralLargestEig.h:445
mat::arn::LanczosSeveralLargestEig::unset_use_full_orth
void unset_use_full_orth()
Definition: LanczosSeveralLargestEig.h:106
mat::arn::LanczosSeveralLargestEig::unset_use_selective_orth
void unset_use_selective_orth()
Definition: LanczosSeveralLargestEig.h:105
mat::arn::LanczosSeveralLargestEig::converged_ith
virtual bool converged_ith(int i) const
Definition: LanczosSeveralLargestEig.h:431
mat::transpose
Xtrans< TX > transpose(TX const &A)
Transposition.
Definition: matrix_proxy.h:131
mat::arn::LanczosSeveralLargestEig::eigVectorTri
Treal * eigVectorTri
Definition: LanczosSeveralLargestEig.h:187
mat
Definition: allocate.cc:39
A
#define A
mat::arn::LanczosSeveralLargestEig::relTol
Treal relTol
Definition: LanczosSeveralLargestEig.h:194
mat::arn::LanczosSeveralLargestEig::size_accTmp
int size_accTmp
Definition: LanczosSeveralLargestEig.h:206
mat::arn::LanczosSeveralLargestEig::setAbsTol
void setAbsTol(Treal const newTol)
Definition: LanczosSeveralLargestEig.h:101
mat::arn::LanczosSeveralLargestEig::counter_orth
int counter_orth
Definition: LanczosSeveralLargestEig.h:217
mat::arn::LanczosSeveralLargestEig::get_num_iter
int get_num_iter() const
Definition: LanczosSeveralLargestEig.h:159
max
#define max(a, b)
Definition: integrator.cc:87
mat::arn::LanczosSeveralLargestEig::run
virtual void run()
Definition: LanczosSeveralLargestEig.h:108
mat::arn::LanczosSeveralLargestEig::~LanczosSeveralLargestEig
virtual ~LanczosSeveralLargestEig()
Definition: LanczosSeveralLargestEig.h:161
mat::arn::LanczosSeveralLargestEig::use_selective_orth
bool use_selective_orth
Definition: LanczosSeveralLargestEig.h:213
mat::arn::LanczosSeveralLargestEig::step
virtual void step()
Definition: LanczosSeveralLargestEig.h:271
mat::arn::LanczosSeveralLargestEig::accTmp
Treal * accTmp
Definition: LanczosSeveralLargestEig.h:204
mat::arn::LanczosSeveralLargestEig::update
virtual void update()
Definition: LanczosSeveralLargestEig.h:197
mat::arn::LanczosSeveralLargestEig::r
Tvector r
Vectors spanning Krylov subspace.
Definition: LanczosSeveralLargestEig.h:185
mat::arn::LanczosSeveralLargestEig::selective_orth
void selective_orth()
Definition: LanczosSeveralLargestEig.h:226
mat::arn::LanczosSeveralLargestEig::getEigVector
void getEigVector(Tvector &eigVec, Treal const *const eVecTri) const
Definition: LanczosSeveralLargestEig.h:390
mat::arn::LanczosSeveralLargestEig::alpha
Treal alpha
Definition: LanczosSeveralLargestEig.h:208
mat::arn::LanczosSeveralLargestEig::eValTmp
Treal * eValTmp
Definition: LanczosSeveralLargestEig.h:203
mat::arn::LanczosSeveralLargestEig::v
Tvector * v
Definition: LanczosSeveralLargestEig.h:180
mat::arn::LanczosSeveralLargestEig::get_ith_eigenpair
virtual void get_ith_eigenpair(int i, Treal &eigVal, Tvector &eigVec, Treal &acc)
Definition: LanczosSeveralLargestEig.h:149
mat::arn::LanczosSeveralLargestEig::set_use_full_orth
void set_use_full_orth()
Definition: LanczosSeveralLargestEig.h:104