ergo
LanczosLargestMagnitudeEig.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 
39 #ifndef MAT_LANCZOSLARGESTMAGNITUDEEIG
40 #define MAT_LANCZOSLARGESTMAGNITUDEEIG
41 #include <limits>
42 #include "Lanczos.h"
43 namespace 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() {
68  computeEigVec();
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;
102  computeEigVec();
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;
117  Treal* eigVectorTri;
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>
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)
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;
227  virtual void update() {
230  }
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
static Treal getRelPrecision()
Definition: matInclude.h:152
Treal eValTmp
Definition: LanczosLargestMagnitudeEig.h:128
void getLargestMagnitudeEigPair(Treal &eValue, Tvector &eVector, Treal &accuracy)
Definition: LanczosLargestMagnitudeEig.h:135
void setAbsTol(Treal const newTol)
Definition: LanczosLargestMagnitudeEig.h:56
Treal absTol
Eigenvector to the tridiagonal matrix length: this->j.
Definition: LanczosLargestMagnitudeEig.h:120
virtual void run()
Definition: Lanczos.h:87
Treal * eigVectorTri
Definition: LanczosLargestMagnitudeEig.h:117
mat::SizesAndBlocks cols
Definition: test.cc:52
Definition: LanczosLargestMagnitudeEig.h:203
Class for building Krylov subspaces with the Lanczos method.
Definition: LanczosLargestMagnitudeEig.h:46
Tvector eVec
Definition: LanczosLargestMagnitudeEig.h:115
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
Treal acc
Definition: LanczosLargestMagnitudeEig.h:116
virtual bool converged() const
Definition: LanczosLargestMagnitudeEig.h:237
virtual void update()
Definition: LanczosLargestMagnitudeEig.h:227
bool eigIsSmall
Definition: LanczosLargestMagnitudeEig.h:226
Definition: allocate.cc:39
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:45
void computeEigenPairTri()
Definition: LanczosLargestMagnitudeEig.h:146
generalVector VectorType
Definition: GetDensFromFock.cc:62
LanczosLargestMagnitudeEigIfSmall(Tmatrix const &AA, Tvector const &startVec, Treal const maxAbsVal, int maxIter=100, int cap=100)
Definition: LanczosLargestMagnitudeEig.h:207
Treal template_blas_fabs(Treal x)
Definition: Interval.h:46
void computeEigVec()
Definition: LanczosLargestMagnitudeEig.h:180
LanczosLargestMagnitudeEig(Tmatrix const &AA, Tvector const &startVec, int maxIter=100, int cap=100)
Definition: LanczosLargestMagnitudeEig.h:49
virtual bool converged() const
Definition: LanczosLargestMagnitudeEig.h:189
void rerun()
Definition: LanczosLargestMagnitudeEig.h:84
Treal eVal
Definition: LanczosLargestMagnitudeEig.h:114
Treal accTmp
Definition: LanczosLargestMagnitudeEig.h:129
int maxIter
Current step.
Definition: Lanczos.h:114
bool largestMagEigIsSmall()
Definition: LanczosLargestMagnitudeEig.h:213
static Treal template_blas_get_num_limit_max()
Definition: template_blas_num_limits.h:85
virtual void run()
Definition: LanczosLargestMagnitudeEig.h:215
Tmatrix const & A
Definition: Lanczos.h:104
void getLargestMagnitudeEig(Treal &ev, Treal &accuracy)
Definition: LanczosLargestMagnitudeEig.h:58
virtual ~LanczosLargestMagnitudeEig()
Definition: LanczosLargestMagnitudeEig.h:110
void restart(Tvector const &startVec)
Definition: Lanczos.h:74
Treal relTol
Definition: LanczosLargestMagnitudeEig.h:121
Class template for building Krylov subspaces with Lanczos.
Definition: Lanczos.h:59
Treal const maxAbsValue
Definition: LanczosLargestMagnitudeEig.h:225
virtual void run()
Definition: LanczosLargestMagnitudeEig.h:66
Xtrans< TX > transpose(TX const &A)
Transposition.
Definition: matrix_proxy.h:131
void setRelTol(Treal const newTol)
Definition: LanczosLargestMagnitudeEig.h:55
virtual void update()
Definition: LanczosLargestMagnitudeEig.h:124
Treal template_blas_sqrt(Treal x)