ergo
LanczosLargestMagnitudeEig.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 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
39#ifndef MAT_LANCZOSLARGESTMAGNITUDEEIG
40#define MAT_LANCZOSLARGESTMAGNITUDEEIG
41#include <limits>
42#include "Lanczos.h"
43namespace mat { /* Matrix namespace */
44 namespace arn { /* Arnoldi type methods namespace */
45 template<typename Treal, typename Tmatrix, typename Tvector>
47 : public Lanczos<Treal, Tmatrix, Tvector> {
48 public:
49 LanczosLargestMagnitudeEig(Tmatrix const & AA, Tvector const & startVec,
50 int maxIter = 100, int cap = 100)
51 : Lanczos<Treal, Tmatrix, Tvector>(AA, startVec, maxIter, cap),
54 eValTmp(0), accTmp(0) {}
55 void setRelTol(Treal const newTol) { relTol = newTol; }
56 void setAbsTol(Treal const newTol) { absTol = newTol; }
57
58 inline void getLargestMagnitudeEig(Treal& ev, Treal& accuracy) {
59 ev = eVal;
60 accuracy = acc;
61 }
62 void getLargestMagnitudeEigPair(Treal& eValue,
63 Tvector& eVector,
64 Treal& accuracy);
65
66 virtual void run() {
69 eVal = eValTmp;
70 acc = accTmp;
71 rerun();
72 /* FIXME! Elias added these extra Lanczos reruns for small
73 matrices to make the tests pass. This is bad.
74 Elias note 2011-01-19: now added one more rerun()
75 because otherwise the test failed when run on Mac.
76 This is really bad. */
77 if(this->A.get_nrows() == 5) {
78 rerun();
79 rerun();
80 rerun();
81 }
82 }
83
84 void rerun() {
85#if 1
86 /* Because of the statistical chance of misconvergence:
87 * Compute new eigenpair with eigenvector in direction orthogonal
88 * to the first eigenvector.
89 */
90 Tvector newResidual(eVec);
91 newResidual.rand();
92 Treal sP = transpose(eVec) * newResidual;
93 newResidual += -sP * eVec;
94 this->restart(newResidual);
96 /* If the new eigenvalue has larger magnitude:
97 * Use those instead
98 */
100 eVal = eValTmp;
101 acc = accTmp;
103 }
104 // Guard against unrealistically small accuracies
105 Treal machine_epsilon = mat::getMachineEpsilon<Treal>();
106 acc = acc / template_blas_fabs(eVal) > 2 * machine_epsilon ? acc : 2 * template_blas_fabs(eVal) * machine_epsilon;
107#endif
108 }
109
111 delete[] eigVectorTri;
112 }
113 protected:
114 Treal eVal;
115 Tvector eVec;
116 Treal acc;
120 Treal absTol;
121 Treal relTol;
122 void computeEigenPairTri();
123 void computeEigVec();
124 virtual void update() {
126 }
127 virtual bool converged() const;
128 Treal eValTmp;
129 Treal accTmp;
130 private:
131 };
132
133 template<typename Treal, typename Tmatrix, typename Tvector>
135 getLargestMagnitudeEigPair(Treal& eValue,
136 Tvector& eVector,
137 Treal& accuracy) {
138 eValue = eVal;
139 accuracy = acc;
140 eVector = eVec;
141 }
142
143
144 template<typename Treal, typename Tmatrix, typename Tvector>
147 delete[] eigVectorTri;
148 Treal* eigVectorMax = new Treal[this->j];
149 Treal* eigVectorMin = new Treal[this->j];
150
151 /* Get largest eigenvalue. */
152 int const lMax = this->j - 1;
153 Treal eValMax;
154 Treal accMax;
155 this->Tri.getEigsByIndex(&eValMax, eigVectorMax, &accMax,
156 lMax, lMax);
157 /* Get smallest eigenvalue. */
158 int const lMin = 0;
159 Treal eValMin;
160 Treal accMin;
161 this->Tri.getEigsByIndex(&eValMin, eigVectorMin, &accMin,
162 lMin, lMin);
163 if (template_blas_fabs(eValMin) > template_blas_fabs(eValMax)) {
164 eValTmp = eValMin;
165 accTmp = accMin;
166 delete[] eigVectorMax;
167 eigVectorTri = eigVectorMin;
168 }
169 else {
170 eValTmp = eValMax;
171 accTmp = accMax;
172 delete[] eigVectorMin;
173 eigVectorTri = eigVectorMax;
174 }
175 }
176
177
178 template<typename Treal, typename Tmatrix, typename Tvector>
181 /* Compute eigenvector as a linear combination of Krylov vectors */
182 assert(eigVectorTri);
183 this->getEigVector(eVec, eigVectorTri);
184 }
185
186
187 template<typename Treal, typename Tmatrix, typename Tvector>
189 converged() const {
190 bool conv = accTmp < absTol; /* Absolute accuracy */
191 if (template_blas_fabs(eValTmp) > 0) {
192 conv = conv &&
193 accTmp / template_blas_fabs(eValTmp) < relTol; /* Relative acc.*/
194 }
195 return conv;
196 }
197
198
199
200
201#if 1
202 template<typename Treal, typename Tmatrix, typename Tvector>
204 : public LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector> {
205 public:
207 (Tmatrix const & AA, Tvector const & startVec,
208 Treal const maxAbsVal,
209 int maxIter = 100, int cap = 100)
210 : LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector>
211 (AA, startVec, maxIter, cap), maxAbsValue(maxAbsVal) {
212 }
214
215 virtual void run() {
217 this->computeEigVec();
218 this->eVal = this->eValTmp;
219 this->acc = this->accTmp;
220 if (eigIsSmall) /* No need to rerun if eigenvalue is large. */
221 this->rerun();
222 }
223
224 protected:
225 Treal const maxAbsValue;
231 virtual bool converged() const;
232 private:
233 };
234
235 template<typename Treal, typename Tmatrix, typename Tvector>
237 converged() const {
238 bool convAccuracy =
240 return convAccuracy || (!eigIsSmall);
241 }
242
243#endif
244
245 } /* end namespace arn */
246
248
249 // EMANUEL COMMENT:
250 // A function similar to euclIfSmall below lies/used to lie inside the MatrixSymmetric class.
251 // There, the matrix was copied and truncated before the calculation.
252 // It is unclear if that had a significant positive impact on the execution time - it definitely required more memory.
259 template<typename Tmatrix, typename Treal>
260 Interval<Treal> euclIfSmall(Tmatrix const & M,
261 Treal const requestedAbsAccuracy,
262 Treal const requestedRelAccuracy,
263 Treal const maxAbsVal,
264 typename Tmatrix::VectorType * const eVecPtr = 0,
265 int maxIter = -1) {
266 assert(requestedAbsAccuracy >= 0);
267 // Treal frobTmp;
268 /* Check if norm is really small, or can easily be determined to be 'large', in those cases quick return */
269 Treal euclLowerBound;
270 Treal euclUpperBound;
271 M.quickEuclBounds(euclLowerBound, euclUpperBound);
272 if ( euclUpperBound < requestedAbsAccuracy )
273 // Norm is really small, quick return
274 return Interval<Treal>( euclLowerBound, euclUpperBound );
275 if ( euclLowerBound > maxAbsVal )
276 // Norm is not small, quick return
277 return Interval<Treal>( euclLowerBound, euclUpperBound );
278 if(maxIter == -1)
279 maxIter = M.get_nrows() * 100;
280 typename Tmatrix::VectorType guess;
282 M.getCols(cols);
283 guess.resetSizesAndBlocks(cols);
284 guess.rand();
286 lan(M, guess, maxAbsVal, maxIter);
287 lan.setAbsTol( requestedAbsAccuracy );
288 lan.setRelTol( requestedRelAccuracy );
289 lan.run();
290 Treal eVal = 0;
291 Treal acc = 0;
292 Treal eValMin = 0;
293 if (lan.largestMagEigIsSmall()) {
294 if (eVecPtr)
295 lan.getLargestMagnitudeEigPair(eVal, *eVecPtr, acc);
296 else
297 lan.getLargestMagnitudeEig(eVal, acc);
298 eVal = template_blas_fabs(eVal);
299 eValMin = eVal - acc;
300 eValMin = eValMin >= 0 ? eValMin : (Treal)0;
301 return Interval<Treal>(eValMin, eVal + acc);
302 }
303 else {
304 eValMin = euclLowerBound;
305 eValMin = eValMin >= maxAbsVal ? eValMin : maxAbsVal;
306 return Interval<Treal>(eValMin, euclUpperBound);
307 }
308 }
309
310} /* end namespace mat */
311#endif
Class for building Krylov subspaces with the Lanczos method.
Definition Interval.h:46
Describes dimensions of matrix and its blocks on all levels.
Definition SizesAndBlocks.h:45
Definition LanczosLargestMagnitudeEig.h:204
virtual void run()
Definition LanczosLargestMagnitudeEig.h:215
virtual bool converged() const
Definition LanczosLargestMagnitudeEig.h:237
LanczosLargestMagnitudeEigIfSmall(Tmatrix const &AA, Tvector const &startVec, Treal const maxAbsVal, int maxIter=100, int cap=100)
Definition LanczosLargestMagnitudeEig.h:207
bool largestMagEigIsSmall()
Definition LanczosLargestMagnitudeEig.h:213
virtual void update()
Definition LanczosLargestMagnitudeEig.h:227
bool eigIsSmall
Definition LanczosLargestMagnitudeEig.h:226
Treal const maxAbsValue
Definition LanczosLargestMagnitudeEig.h:225
Definition LanczosLargestMagnitudeEig.h:47
void computeEigenPairTri()
Definition LanczosLargestMagnitudeEig.h:146
void setRelTol(Treal const newTol)
Definition LanczosLargestMagnitudeEig.h:55
Treal eValTmp
Definition LanczosLargestMagnitudeEig.h:128
void setAbsTol(Treal const newTol)
Definition LanczosLargestMagnitudeEig.h:56
Tvector eVec
Definition LanczosLargestMagnitudeEig.h:115
virtual void run()
Definition LanczosLargestMagnitudeEig.h:66
Treal eVal
Definition LanczosLargestMagnitudeEig.h:114
virtual ~LanczosLargestMagnitudeEig()
Definition LanczosLargestMagnitudeEig.h:110
LanczosLargestMagnitudeEig(Tmatrix const &AA, Tvector const &startVec, int maxIter=100, int cap=100)
Definition LanczosLargestMagnitudeEig.h:49
void getLargestMagnitudeEigPair(Treal &eValue, Tvector &eVector, Treal &accuracy)
Definition LanczosLargestMagnitudeEig.h:135
Treal absTol
Eigenvector to the tridiagonal matrix length: this->j.
Definition LanczosLargestMagnitudeEig.h:120
Treal acc
Definition LanczosLargestMagnitudeEig.h:116
void getLargestMagnitudeEig(Treal &ev, Treal &accuracy)
Definition LanczosLargestMagnitudeEig.h:58
virtual bool converged() const
Definition LanczosLargestMagnitudeEig.h:189
virtual void update()
Definition LanczosLargestMagnitudeEig.h:124
Treal relTol
Definition LanczosLargestMagnitudeEig.h:121
void rerun()
Definition LanczosLargestMagnitudeEig.h:84
void computeEigVec()
Definition LanczosLargestMagnitudeEig.h:180
Treal accTmp
Definition LanczosLargestMagnitudeEig.h:129
Treal * eigVectorTri
Definition LanczosLargestMagnitudeEig.h:117
Class template for building Krylov subspaces with Lanczos.
Definition Lanczos.h:59
virtual void run()
Definition Lanczos.h:87
void restart(Tvector const &startVec)
Definition Lanczos.h:74
int maxIter
Current step.
Definition Lanczos.h:114
Tmatrix const & A
Definition Lanczos.h:104
mat::SizesAndBlocks cols
Definition test.cc:52
Definition allocate.cc:39
static Treal getMachineEpsilon()
Definition matInclude.h:147
Xtrans< TX > transpose(TX const &A)
Transposition.
Definition matrix_proxy.h:131
Interval< Treal > euclIfSmall(Tmatrix const &M, Treal const requestedAbsAccuracy, Treal const requestedRelAccuracy, Treal const maxAbsVal, typename Tmatrix::VectorType *const eVecPtr=0, int maxIter=-1)
Returns interval containing the Euclidean norm of the matrix M.
Definition LanczosLargestMagnitudeEig.h:260
static Treal getRelPrecision()
Definition matInclude.h:152
Treal template_blas_sqrt(Treal x)
Treal template_blas_fabs(Treal x)
static Treal template_blas_get_num_limit_max()
Definition template_blas_num_limits.h:85