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
46namespace 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)
112 step();
113 update();
114 if (j > maxIter)
115 throw AcceptableMaxIter("Lanczos::run() did not converge within maxIter");
116 }
117 while (!converged());
118
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
189 int j;
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
Definition: Failure.h:81
Definition: LanczosSeveralLargestEig.h:51
Treal relTol
Definition: LanczosSeveralLargestEig.h:194
Tvector * deflVec
Definition: LanczosSeveralLargestEig.h:220
void unset_use_full_orth()
Definition: LanczosSeveralLargestEig.h:106
void selective_orth()
Definition: LanczosSeveralLargestEig.h:226
MatrixTridiagSymmetric< Treal > Tri
Residual vector.
Definition: LanczosSeveralLargestEig.h:186
void copyTridiag(MatrixTridiagSymmetric< Treal > &Tricopy)
Definition: LanczosSeveralLargestEig.h:173
int number_of_eigenv
Definition: LanczosSeveralLargestEig.h:205
virtual bool converged_ith(int i) const
Definition: LanczosSeveralLargestEig.h:431
Treal * accTmp
Definition: LanczosSeveralLargestEig.h:204
Treal beta
Definition: LanczosSeveralLargestEig.h:209
Treal * eValTmp
Definition: LanczosSeveralLargestEig.h:203
int get_num_iter() const
Definition: LanczosSeveralLargestEig.h:159
bool use_selective_orth
Definition: LanczosSeveralLargestEig.h:213
int counter_orth
Definition: LanczosSeveralLargestEig.h:217
Treal * eigVectorTri
Definition: LanczosSeveralLargestEig.h:187
virtual bool converged() const
Definition: LanczosSeveralLargestEig.h:417
void unset_use_selective_orth()
Definition: LanczosSeveralLargestEig.h:105
Tmatrix const & A
Definition: LanczosSeveralLargestEig.h:179
virtual ~LanczosSeveralLargestEig()
Definition: LanczosSeveralLargestEig.h:161
virtual void step()
Definition: LanczosSeveralLargestEig.h:271
virtual void run()
Definition: LanczosSeveralLargestEig.h:108
void setAbsTol(Treal const newTol)
Definition: LanczosSeveralLargestEig.h:101
void set_use_full_orth()
Definition: LanczosSeveralLargestEig.h:104
void set_use_selective_orth()
Definition: LanczosSeveralLargestEig.h:103
Tvector r
Vectors spanning Krylov subspace.
Definition: LanczosSeveralLargestEig.h:185
void getEigVector(Tvector &eigVec, Treal const *const eVecTri) const
Definition: LanczosSeveralLargestEig.h:390
Treal alpha
Definition: LanczosSeveralLargestEig.h:208
Treal sigma
Definition: LanczosSeveralLargestEig.h:221
Treal absTol
Definition: LanczosSeveralLargestEig.h:193
virtual void computeEigenPairTri()
Definition: LanczosSeveralLargestEig.h:343
bool use_full_orth
Definition: LanczosSeveralLargestEig.h:214
virtual void get_ith_eigenpair(int i, Treal &eigVal, Tvector &eigVec, Treal &acc)
Definition: LanczosSeveralLargestEig.h:149
Tvector * v
Definition: LanczosSeveralLargestEig.h:180
int capacity
Definition: LanczosSeveralLargestEig.h:188
void setRelTol(Treal const newTol)
Definition: LanczosSeveralLargestEig.h:100
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 total_num_iter
Definition: LanczosSeveralLargestEig.h:211
int counter_all
Definition: LanczosSeveralLargestEig.h:216
void increaseCapacity(int const newCapacity)
Definition: LanczosSeveralLargestEig.h:445
int size_accTmp
Definition: LanczosSeveralLargestEig.h:206
virtual void update()
Definition: LanczosSeveralLargestEig.h:197
int maxIter
Current step.
Definition: LanczosSeveralLargestEig.h:190
int j
Definition: LanczosSeveralLargestEig.h:189
Tridiagonal symmetric matrix class template.
Definition: MatrixTridiagSymmetric.h:47
#define max(a, b)
Definition: integrator.cc:87
#define A
Definition: allocate.cc:39
Xtrans< TX > transpose(TX const &A)
Transposition.
Definition: matrix_proxy.h:131
Treal template_blas_sqrt(Treal x)
Treal template_blas_fabs(Treal x)