ergo
LanczosLargestMagnitudeEig.h
Go to the documentation of this file.
1 /* Ergo, version 3.2, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
37 #ifndef MAT_LANCZOSLARGESTMAGNITUDEEIG
38 #define MAT_LANCZOSLARGESTMAGNITUDEEIG
39 #include <limits>
40 #include "Lanczos.h"
41 namespace mat { /* Matrix namespace */
42  namespace arn { /* Arnoldi type methods namespace */
43  template<typename Treal, typename Tmatrix, typename Tvector>
45  : public Lanczos<Treal, Tmatrix, Tvector> {
46  public:
47  LanczosLargestMagnitudeEig(Tmatrix const & AA, Tvector const & startVec,
48  int maxIter = 100, int cap = 100)
49  : Lanczos<Treal, Tmatrix, Tvector>(AA, startVec, maxIter, cap),
50  eVal(0), acc(0), eigVectorTri(0), absTol( std::numeric_limits<Treal>::max() ),
52  eValTmp(0), accTmp(0) {}
53  void setRelTol(Treal const newTol) { relTol = newTol; }
54  void setAbsTol(Treal const newTol) { absTol = newTol; }
55 
56  inline void getLargestMagnitudeEig(Treal& ev, Treal& accuracy) {
57  ev = eVal;
58  accuracy = acc;
59  }
60  void getLargestMagnitudeEigPair(Treal& eValue,
61  Tvector& eVector,
62  Treal& accuracy);
63 
64  virtual void run() {
66  computeEigVec();
67  eVal = eValTmp;
68  acc = accTmp;
69  rerun();
70  /* FIXME! Elias added these extra Lanczos reruns for small
71  matrices to make the tests pass. This is bad.
72  Elias note 2011-01-19: now added one more rerun()
73  because otherwise the test failed when run on Mac.
74  This is really bad. */
75  if(this->A.get_nrows() == 5) {
76  rerun();
77  rerun();
78  rerun();
79  }
80  }
81 
82  void rerun() {
83 #if 1
84  /* Because of the statistical chance of misconvergence:
85  * Compute new eigenpair with eigenvector in direction orthogonal
86  * to the first eigenvector.
87  */
88  Tvector newResidual(eVec);
89  newResidual.rand();
90  Treal sP = transpose(eVec) * newResidual;
91  newResidual += -sP * eVec;
92  this->restart(newResidual);
94  /* If the new eigenvalue has larger magnitude:
95  * Use those instead
96  */
98  eVal = eValTmp;
99  acc = accTmp;
100  computeEigVec();
101  }
102  // Guard against unrealistically small accuracies
103  acc = acc / template_blas_fabs(eVal) > 2 * std::numeric_limits<Treal>::epsilon() ? acc : 2 * template_blas_fabs(eVal) * std::numeric_limits<Treal>::epsilon();
104 #endif
105  }
106 
108  delete[] eigVectorTri;
109  }
110  protected:
111  Treal eVal;
112  Tvector eVec;
113  Treal acc;
114  Treal* eigVectorTri;
117  Treal absTol;
118  Treal relTol;
119  void computeEigenPairTri();
120  void computeEigVec();
121  virtual void update() {
123  }
124  virtual bool converged() const;
125  Treal eValTmp;
126  Treal accTmp;
127  private:
128  };
129 
130  template<typename Treal, typename Tmatrix, typename Tvector>
133  Tvector& eVector,
134  Treal& accuracy) {
135  eValue = eVal;
136  accuracy = acc;
137  eVector = eVec;
138  }
139 
140 
141  template<typename Treal, typename Tmatrix, typename Tvector>
144  delete[] eigVectorTri;
145  Treal* eigVectorMax = new Treal[this->j];
146  Treal* eigVectorMin = new Treal[this->j];
147 
148  /* Get largest eigenvalue. */
149  int const lMax = this->j - 1;
150  Treal eValMax;
151  Treal accMax;
152  this->Tri.getEigsByIndex(&eValMax, eigVectorMax, &accMax,
153  lMax, lMax);
154  /* Get smallest eigenvalue. */
155  int const lMin = 0;
156  Treal eValMin;
157  Treal accMin;
158  this->Tri.getEigsByIndex(&eValMin, eigVectorMin, &accMin,
159  lMin, lMin);
160  if (template_blas_fabs(eValMin) > template_blas_fabs(eValMax)) {
161  eValTmp = eValMin;
162  accTmp = accMin;
163  delete[] eigVectorMax;
164  eigVectorTri = eigVectorMin;
165  }
166  else {
167  eValTmp = eValMax;
168  accTmp = accMax;
169  delete[] eigVectorMin;
170  eigVectorTri = eigVectorMax;
171  }
172  }
173 
174 
175  template<typename Treal, typename Tmatrix, typename Tvector>
178  /* Compute eigenvector as a linear combination of Krylov vectors */
179  assert(eigVectorTri);
180  this->getEigVector(eVec, eigVectorTri);
181  }
182 
183 
184  template<typename Treal, typename Tmatrix, typename Tvector>
186  converged() const {
187  bool conv = accTmp < absTol; /* Absolute accuracy */
188  if (template_blas_fabs(eValTmp) > 0) {
189  conv = conv &&
190  accTmp / template_blas_fabs(eValTmp) < relTol; /* Relative acc.*/
191  }
192  return conv;
193  }
194 
195 
196 
197 
198 #if 1
199  template<typename Treal, typename Tmatrix, typename Tvector>
201  : public LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector> {
202  public:
204  (Tmatrix const & AA, Tvector const & startVec,
205  Treal const maxAbsVal,
206  int maxIter = 100, int cap = 100)
208  (AA, startVec, maxIter, cap), maxAbsValue(maxAbsVal) {
209  }
211 
212  virtual void run() {
214  this->computeEigVec();
215  this->eVal = this->eValTmp;
216  this->acc = this->accTmp;
217  if (eigIsSmall) /* No need to rerun if eigenvalue is large. */
218  this->rerun();
219  }
220 
221  protected:
222  Treal const maxAbsValue;
224  virtual void update() {
227  }
228  virtual bool converged() const;
229  private:
230  };
231 
232  template<typename Treal, typename Tmatrix, typename Tvector>
234  converged() const {
235  bool convAccuracy =
237  return convAccuracy || (!eigIsSmall);
238  }
239 
240 #endif
241 
242  } /* end namespace arn */
243 
245 
246  // EMANUEL COMMENT:
247  // A function similar to euclIfSmall below lies/used to lie inside the MatrixSymmetric class.
248  // There, the matrix was copied and truncated before the calculation.
249  // It is unclear if that had a significant positive impact on the execution time - it definitely required more memory.
256  template<typename Tmatrix, typename Treal>
257  Interval<Treal> euclIfSmall(Tmatrix const & M,
258  Treal const requestedAbsAccuracy,
259  Treal const requestedRelAccuracy,
260  Treal const maxAbsVal,
261  typename Tmatrix::VectorType * const eVecPtr = 0) {
262  assert(requestedAbsAccuracy >= 0);
263  // Treal frobTmp;
264  /* Check if norm is really small, or can easily be determined to be 'large', in those cases quick return */
265  Treal euclLowerBound;
266  Treal euclUpperBound;
267  M.quickEuclBounds(euclLowerBound, euclUpperBound);
268  if ( euclUpperBound < requestedAbsAccuracy )
269  // Norm is really small, quick return
270  return Interval<Treal>( euclLowerBound, euclUpperBound );
271  if ( euclLowerBound > maxAbsVal )
272  // Norm is not small, quick return
273  return Interval<Treal>( euclLowerBound, euclUpperBound );
274  int maxIter = M.get_nrows() * 100;
275  typename Tmatrix::VectorType guess;
276  SizesAndBlocks cols;
277  M.getCols(cols);
278  guess.resetSizesAndBlocks(cols);
279  guess.rand();
281  lan(M, guess, maxAbsVal, maxIter);
282  lan.setAbsTol( requestedAbsAccuracy );
283  lan.setRelTol( requestedRelAccuracy );
284  lan.run();
285  Treal eVal = 0;
286  Treal acc = 0;
287  Treal eValMin = 0;
288  if (lan.largestMagEigIsSmall()) {
289  if (eVecPtr)
290  lan.getLargestMagnitudeEigPair(eVal, *eVecPtr, acc);
291  else
292  lan.getLargestMagnitudeEig(eVal, acc);
293  eVal = template_blas_fabs(eVal);
294  eValMin = eVal - acc;
295  eValMin = eValMin >= 0 ? eValMin : (Treal)0;
296  return Interval<Treal>(eValMin, eVal + acc);
297  }
298  else {
299  eValMin = euclLowerBound;
300  eValMin = eValMin >= maxAbsVal ? eValMin : maxAbsVal;
301  return Interval<Treal>(eValMin, euclUpperBound);
302  }
303  }
304 
305 } /* end namespace mat */
306 #endif