ergo
LanczosSeveralLargestEig.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 
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 
120 
121 
122 
123  // i is a number of eigenvalue (1 is the largest, 2 is the second largest and so on)
124  virtual void get_ith_eigenpair(int i, Treal& eigVal, Tvector& eigVec, Treal & acc)
125  {
126  assert(i > 0);
127  assert(i <= size_accTmp);
128  eigVal = eValTmp[size_accTmp - i]; // array
129  assert(eigVectorTri);
130  getEigVector(eigVec, &eigVectorTri[j * (size_accTmp - i)]);
131  acc = accTmp[size_accTmp - i];
132  }
133 
134  int get_num_iter() const{ return j;}
135 
137 
139  printf("Orthogonalized %d of total possible %d, this is %lf %%\n", counter_orth, counter_all, (double)counter_orth/counter_all*100);
140 
141  delete[] eigVectorTri;
142  delete[] eValTmp;
143  delete[] accTmp;
144  delete[] v;
145  }
146 
147 
149  Tricopy = Tri;
150  }
151 
152 
153  protected:
154  Tmatrix const & A;
155  Tvector* v;
160  Tvector r;
162  Treal* eigVectorTri; // Eigenvectors of the tridiagonal matrix
163  int capacity;
164  int j;
165  int maxIter;
166  void increaseCapacity(int const newCapacity);
167  void getEigVector(Tvector& eigVec, Treal const * const eVecTri) const;
168  Treal absTol;
169  Treal relTol;
170  virtual void step();
171  virtual void computeEigenPairTri();
172  virtual void update() {
174  }
175  void selective_orth();
176  virtual bool converged() const;
177  virtual bool converged_ith(int i) const;
178  Treal* eValTmp; // current computed eigenvalues (less or equal to number_of_eigenv)
179  Treal* accTmp; // residuals
180  int number_of_eigenv; // eigenvalues are saved in the decreasing order, thus the largest one has index 1
181  int size_accTmp; // size of accTmp (number of computed eigenvalues of the matrix T)
182  private:
183  Treal alpha;
184  Treal beta;
185 
188 
191 
192  // if deflation is used
193  Tvector * deflVec;
194  Treal sigma;
195  };
196 
197  template<typename Treal, typename Tmatrix, typename Tvector>
200  {
201  int j_curr = j-1;
202 
203  Treal coeff = 0, res;
204  Treal normT = 0; // spectral norm of T (since norm of A is not available)
205  // find largest by absolute value eigenvalue of T
206  for(int i = 0; i <= j_curr; ++i)
207  if(template_blas_fabs(eValTmp[i]) > normT) normT = template_blas_fabs(eValTmp[i]);
208 
209  Treal epsilon = mat::getMachineEpsilon<Treal>();
210  Tvector tmp;
211  tmp = v[j_curr+1];
212  tmp *= beta; // return non-normalized value
213 
214  for(int i = j_curr; i >= 0; --i)
215  {
216  counter_all++;
217  // get residual for this eigenpair
218  res = accTmp[i];
219  Treal tol = template_blas_sqrt(epsilon) * normT;
220  if(res <= tol) // b_{j} * |VT_i(j)| <= sqrt(eps) * norm(A), but we do not have norm(A)
221  {
222  counter_orth++;
223  Tvector eigVec;
224  getEigVector(eigVec, &eigVectorTri[j_curr * i]); // y = U*VT(:, i); % ith Ritz vector
225  coeff = transpose(eigVec) * tmp;
226  tmp += (-coeff) * (eigVec); // v = v - (y'*v)*y
227  }
228  }
229 
230 
231 
232  v[j_curr+1] = tmp;
233  beta = v[j_curr+1].eucl(); // update beta
234  Treal const ONE = 1.0;
235  v[j_curr+1] *= ONE / beta; // normalized
236  Tri.update_beta(beta);
237 
238  /* /\* // check orthogonality *\/ */
239  /* try */
240  /* { */
241  /* for(int k = 0; k < j_curr; ++k) */
242  /* { */
243  /* v[k].readFromFile(); */
244  /* Treal val = transpose(v[k]) * v[j_curr+1]; */
245  /* std::cout << val << ", "; */
246  /* v[k].writeToFile(); */
247  /* } */
248  /* std::cout << std::endl; */
249  /* } */
250  /* catch(const std::exception &e) */
251  /* { */
252  /* std::cout << "Exception: " << e.what() << std::endl; */
253  /* } */
254  }
255 
256 
257 
258 
259  template<typename Treal, typename Tmatrix, typename Tvector>
262  {
263  if (j + 1 >= capacity)
264  increaseCapacity(capacity * 2);
265  Treal const ONE = 1.0;
266  A.matVecProd(r, v[j]); // r = A * v[j]
267  alpha = transpose(v[j]) * r; // alpha = v[j]'*A*v[j]
268 
269  /*
270  If one wants to use deflation with vector
271  x_1:=deflVec (usually it is an eigenvector
272  corresponding to an eigenvalue lambda_1 of A)
273  and thus compute eigenvalues of the matrix
274  An = A-sigma*x_1*x_1'
275  Note: if lambga_i are eigenvalues of A corresponding to x_i, then
276  An will have eigenvalues (lambda_1-sigma, lambda_2, ..., lambda_N)
277  and unchanged eigenvectors x_i.
278  */
279 
280  if(deflVec != NULL)
281  {
282  /*
283  r = (A*vj - sigma*(x_1'*vj)*x_1) - alpha*vj - beta*v{j-1}
284  where
285  alpha = vj'*An*vj = vj'*A*vj - sigma * (x_1'*vj)^2
286  */
287  Treal gamma = transpose(*deflVec) * v[j]; // dot product x' * v_j
288  alpha -= sigma*gamma*gamma;
289 
290  r += (-sigma*gamma) * (*deflVec);
291  }
292 
293  r += (-alpha) * v[j];
294 
295  if (j) {
296  r += (-beta) * v[j-1];
297  v[j-1].writeToFile();
298  }
299 
300 
301 
302  if(use_full_orth)
303  {
304  // full orthogonalization
305  // r = r - (q_1'*r)*q_1 - ... - (q_{j-1}'*r)*q_{j-1}
306  Treal gamma_i = 0;
307  Tvector tmp;
308  tmp = r;
309  for(int i = 0; i < j; ++i )
310  {
311  v[i].readFromFile();
312  gamma_i = transpose(tmp) * v[i]; // r'*v_i
313  r += (-gamma_i) * v[i]; // (r'*vi) * v_i
314  v[i].writeToFile();
315  }
316  tmp.clear();
317  }
318 
319 
320  beta = r.eucl();
321  v[j+1] = r;
322  v[j+1] *= ONE / beta;
323  Tri.increase(alpha, beta);
324 
325 
326  /* /\* // check orthogonality *\/ */
327  /* try */
328  /* { */
329  /* for(int k = 0; k < j; ++k) */
330  /* { */
331  /* v[k].readFromFile(); */
332  /* Treal val = transpose(v[k]) * v[j+1]; */
333  /* std::cout << val << ", "; */
334  /* v[k].writeToFile(); */
335  /* } */
336  /* std::cout << std::endl << "-----" << std::endl; */
337  /* } */
338  /* catch(const std::exception &e) */
339  /* { */
340  /* std::cout << "Exception: " << e.what() << std::endl; */
341  /* } */
342 
343 
344  j++;
345  }
346 
347 
348  /*
349  Compute eigenvectors of the tridiagonal matrix
350  */
351  template<typename Treal, typename Tmatrix, typename Tvector>
354  if( eigVectorTri != NULL ) delete[] eigVectorTri;
355  if( accTmp != NULL ) delete[] accTmp;
356  if( eValTmp != NULL ) delete[] eValTmp;
357 
358  int num_compute_eigenvalues;
359  if(use_selective_orth)
360  num_compute_eigenvalues = j; // we need all eigenvectors of T
361  else
362  num_compute_eigenvalues = number_of_eigenv; // it is enough just number_of_eigenv of T
363 
364  /* Get largest eigenvalues */
365  int const max_ind = j-1; // eigenvalue count starts with 0
366  int const min_ind = std::max(j - num_compute_eigenvalues, 0);
367 
368  Treal* eigVectors = new Treal[j * num_compute_eigenvalues]; // every vector of size j
369  Treal* eigVals = new Treal[num_compute_eigenvalues];
370  Treal* accMax = new Treal[num_compute_eigenvalues];
371  assert(eigVectors != NULL);
372  assert(eigVals != NULL);
373  assert(accMax != NULL);
374 
375  Tri.getEigsByIndex(eigVals, eigVectors, accMax,
376  min_ind, max_ind);
377 
378  eValTmp = eigVals;
379 
380 
381  eigVectorTri = eigVectors;
382  accTmp = accMax;
383  size_accTmp = num_compute_eigenvalues;
384 
385  // set unused pointers to NULL
386  eigVectors = NULL;
387  eigVals = NULL;
388  accMax = NULL;
389  }
390 
391 
392 
393 
394  /* FIXME: If several eigenvectors are needed it is more optimal to
395  * compute all of them at once since then the krylov subspace vectors
396  * only need to be read from memory once.
397  */
398  template<typename Treal, typename Tmatrix, typename Tvector>
400  getEigVector(Tvector& eigVec, Treal const * const eVecTri) const {
401  if (j <= 1) {
402  eigVec = v[0];
403  }
404  else {
405  v[0].readFromFile();
406  eigVec = v[0];
407  v[0].writeToFile();
408  }
409  eigVec *= eVecTri[0];
410  for (int ind = 1; ind <= j - 2; ++ind) {
411  v[ind].readFromFile();
412  eigVec += eVecTri[ind] * v[ind];
413  v[ind].writeToFile();
414  }
415  eigVec += eVecTri[j-1] * v[j-1];
416 
417  // normalized
418  Treal norm_eigVec = eigVec.eucl();
419  Treal const ONE = 1.0;
420  eigVec *= ONE / norm_eigVec;
421  }
422 
423 
424  // we want lowest eigenvalue to converge
425  template<typename Treal, typename Tmatrix, typename Tvector>
427  converged() const {
428 
429  if(j < number_of_eigenv) return false;
430  bool conv = converged_ith(number_of_eigenv); // if the last needed eigenvalue converged
431 
432  return conv;
433  }
434 
435  // check convergence of ith eigenpair
436  template<typename Treal, typename Tmatrix, typename Tvector>
438  converged_ith(int i) const {
439  assert(size_accTmp >= i);
440 
441  bool conv = true; //accTmp[size_accTmp - i] < absTol; /* Do not use absolute accuracy */
442  if (template_blas_fabs(eValTmp[size_accTmp - i]) > 0) {
443  conv = conv &&
444  accTmp[size_accTmp - i] / template_blas_fabs(eValTmp[size_accTmp - i]) < relTol; /* Relative acc.*/
445  }
446  return conv;
447  }
448 
449 
450  template<typename Treal, typename Tmatrix, typename Tvector>
452  increaseCapacity(int const newCapacity) {
453  assert(newCapacity > capacity);
454  assert(j > 0);
455  capacity = newCapacity;
456  Tvector* new_v = new Tvector[capacity];
457  assert(new_v != NULL);
458  /* FIXME: Fix so that file is copied when operator= is called in Vector
459  * class
460  */
461  for (int ind = 0; ind <= j - 2; ind++){
462  v[ind].readFromFile();
463  new_v[ind] = v[ind];
464  new_v[ind].writeToFile();
465  }
466  for (int ind = j - 1; ind <= j; ind++){
467  new_v[ind] = v[ind];
468  }
469  delete[] v;
470  v = new_v;
471  }
472 
473 
474  } /* end namespace arn */
475 
476 
477 } /* end namespace mat */
478 #endif
#define A
Definition: LanczosSeveralLargestEig.h:50
void setRelTol(Treal const newTol)
Definition: LanczosSeveralLargestEig.h:100
Treal sigma
Definition: LanczosSeveralLargestEig.h:194
Tvector * deflVec
Definition: LanczosSeveralLargestEig.h:193
void increaseCapacity(int const newCapacity)
Definition: LanczosSeveralLargestEig.h:452
void unset_use_full_orth()
Definition: LanczosSeveralLargestEig.h:106
Definition: Failure.h:81
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
int counter_all
Definition: LanczosSeveralLargestEig.h:189
Tridiagonal symmetric matrix class template.
Definition: MatrixTridiagSymmetric.h:47
virtual bool converged() const
Definition: LanczosSeveralLargestEig.h:427
virtual void step()
Definition: LanczosSeveralLargestEig.h:261
Tvector * v
Definition: LanczosSeveralLargestEig.h:155
virtual void computeEigenPairTri()
Definition: LanczosSeveralLargestEig.h:353
Definition: allocate.cc:39
int j
Definition: LanczosSeveralLargestEig.h:164
Treal template_blas_fabs(Treal x)
void copyTridiag(MatrixTridiagSymmetric< Treal > &Tricopy)
Definition: LanczosSeveralLargestEig.h:148
Tmatrix const & A
Definition: LanczosSeveralLargestEig.h:154
Treal * eigVectorTri
Definition: LanczosSeveralLargestEig.h:162
#define max(a, b)
Definition: integrator.cc:87
void selective_orth()
Definition: LanczosSeveralLargestEig.h:199
Treal * accTmp
Definition: LanczosSeveralLargestEig.h:179
void set_use_selective_orth()
Definition: LanczosSeveralLargestEig.h:103
Treal * eValTmp
Definition: LanczosSeveralLargestEig.h:178
virtual void get_ith_eigenpair(int i, Treal &eigVal, Tvector &eigVec, Treal &acc)
Definition: LanczosSeveralLargestEig.h:124
virtual ~LanczosSeveralLargestEig()
Definition: LanczosSeveralLargestEig.h:136
void unset_use_selective_orth()
Definition: LanczosSeveralLargestEig.h:105
virtual bool converged_ith(int i) const
Definition: LanczosSeveralLargestEig.h:438
bool use_full_orth
Definition: LanczosSeveralLargestEig.h:187
void getEigVector(Tvector &eigVec, Treal const *const eVecTri) const
Definition: LanczosSeveralLargestEig.h:400
Treal absTol
Definition: LanczosSeveralLargestEig.h:168
Treal relTol
Definition: LanczosSeveralLargestEig.h:169
int size_accTmp
Definition: LanczosSeveralLargestEig.h:181
int get_num_iter() const
Definition: LanczosSeveralLargestEig.h:134
void set_use_full_orth()
Definition: LanczosSeveralLargestEig.h:104
bool use_selective_orth
Definition: LanczosSeveralLargestEig.h:186
void setAbsTol(Treal const newTol)
Definition: LanczosSeveralLargestEig.h:101
MatrixTridiagSymmetric< Treal > Tri
Residual vector.
Definition: LanczosSeveralLargestEig.h:161
Treal alpha
Definition: LanczosSeveralLargestEig.h:183
int maxIter
Current step.
Definition: LanczosSeveralLargestEig.h:165
int number_of_eigenv
Definition: LanczosSeveralLargestEig.h:180
Tvector r
Vectors spanning Krylov subspace.
Definition: LanczosSeveralLargestEig.h:160
virtual void run()
Definition: LanczosSeveralLargestEig.h:108
int counter_orth
Definition: LanczosSeveralLargestEig.h:190
Xtrans< TX > transpose(TX const &A)
Transposition.
Definition: matrix_proxy.h:131
Treal template_blas_sqrt(Treal x)
virtual void update()
Definition: LanczosSeveralLargestEig.h:172
int capacity
Definition: LanczosSeveralLargestEig.h:163
Treal beta
Definition: LanczosSeveralLargestEig.h:184