ergo
Purification.h
Go to the documentation of this file.
1 /* Ergo, version 3.3, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2013 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 
36 #ifndef MAT_PURIFICATION
37 #define MAT_PURIFICATION
38 #include <math.h>
39 #include <iomanip>
40 #define __ID__ "$Id$"
41 namespace mat {
42  template<typename Treal, typename Tmatrix, typename TdebugPolicy>
43  class Purification :public TdebugPolicy {
44  public:
45  typedef typename Tmatrix::VectorType VectorType;
46  Purification(Tmatrix& M,
49  normType const normXmX2,
57  );
58  void step();
59  void purify();
60  protected:
61  Tmatrix & X;
62  Tmatrix X2;
66  int niter;
67 
69 
70  private:
71  };
72 
73  template<typename Treal, typename Tmatrix, typename TdebugPolicy>
75  Purification(Tmatrix& M,
76  normType const normXmX2_inp,
78  :X(M), X2(M), normTruncation( infop.getNormForTruncation() ), normXmX2(normXmX2_inp), info(infop),niter(0) {
79  Time tStepTotal;
80  tStepTotal.tic();
81 
82  Treal lmin(0);
83  Treal lmax(0);
85  currentStep.improveEigInterval(Interval<Treal>(0.0,1.0));
87  lmin = eigFInt.low();
88  lmax = eigFInt.upp();
89  X.add_identity(-lmax); /* Scale to [0, 1] interval and negate */
90  X *= ((Treal)1.0 / (lmin - lmax));
91  /* Compute transformed homo and lumo eigenvalues. */
92  Interval<Treal> homo = info.getHomoF();
93  Interval<Treal> lumo = info.getLumoF();
94  homo = (homo - lmax) / (lmin - lmax);
95  lumo = (lumo - lmax) / (lmin - lmax);
96 
97  Treal chosenThresh = info.getOptimalThresh();
98  currentStep.setChosenThresh(chosenThresh);
99  /* Consider convergence of eigs in getOptimalThresh */
100  currentStep.setMemUsageBeforeTrunc();
101  Time t;
102  t.tic();
103  Treal actualThresh = X.thresh(chosenThresh, normTruncation);
104  currentStep.setTimeThresh(t.toc());
105  currentStep.setActualThresh(actualThresh);
106  if (homo.empty() || lumo.empty()) {
107  throw Failure("Purification<Treal, Tmatrix, TdebugPolicy>::"
108  "Purification(Tmatrix&, normType const, "
109  "PuriInfo<Treal, VectorType, TdebugPolicy>&): "
110  "HOMO or LUMO empty in purification constructor. ");
111  }
112 
113  homo.increase(actualThresh);
114  lumo.increase(actualThresh);
115  stepComputeInfo(currentStep);
116  currentStep.improveHomoLumo(homo,lumo);
117  currentStep.setTimeTotal(tStepTotal.toc());
118  }
120  template<typename Treal, typename Tmatrix, typename TdebugPolicy>
122  Time tStepTotal;
123  tStepTotal.tic();
124  PuriStepInfo<Treal, VectorType, TdebugPolicy> & currentStep = info.getNext();
125  /* It may happen that X2 has many more nonzeros than X, for
126  example 5 times as many. Therefore it makes sense to try
127  having only one "big" matrix in memory at a time. However,
128  file operations have proved to be quite expensive and should
129  be avoided if possible. Hence we want to achieve having only
130  one big matrix in memory without unnecessary file
131  operations. We are currently hoping that it will be ok to add
132  a "small" matrix to a "big" one, that the memory usage after
133  that operation will be like the memory usage for one big
134  matrix + one small matrix. Therefore we are adding X to X2 (X
135  is truncated, a "small" matrix) instead of the opposite when
136  the 2*X-X*X polynomial is evaluated.
137  */
138  if (info(niter).getPoly()) {
139  /* Polynomial 2 * x - x * x */
140  X2 *= (Treal)-1.0;
141  X2 += (Treal)2.0 * X;
142  // Now X2 contains 2*X-X*X
143  }
144  /* In case of polynomial x * x we only need to transfer the
145  content of X2 to X. */
146  // Transfer content of X2 to X clearing previous content of X if any
147  // In current implementation this is needed regardless of which
148  // polynomial is used
149  X2.transfer(X);
150 
151  /* Obtain homo/lumo information from previous before thresh choice.
152  * Default value of 0.0 for thresh is used. getOptimalThresh can
153  * try different thresh values.
154  */
155  Treal chosenThresh = info.getOptimalThresh();
156  currentStep.setChosenThresh(chosenThresh);
157  /* Consider convergence of eigs in getOptimalThresh? */
158  currentStep.setMemUsageBeforeTrunc();
159  Time t;
160  t.tic();
161  currentStep.setActualThresh(X.thresh(chosenThresh, normTruncation));
162  currentStep.setTimeThresh(t.toc());
163  stepComputeInfo(currentStep);
164  if (niter > 5 &&
165  currentStep.getHomo().low() > 0.9 &&
166  currentStep.getLumo().upp() < 0.1 &&
167  info(niter-5).getHomo().low() > 0.9 &&
168  info(niter-5).getLumo().upp() < 0.1 &&
169  ((1 - currentStep.getHomo().low()) >
170  (1 - info(niter-5).getHomo().low()) -
171  currentStep.getEigAccLoss() ||
172  currentStep.getLumo().upp() >
173  info(niter-5).getLumo().upp() -
174  currentStep.getEigAccLoss())) {
175  throw Failure("Purification<Treal, Tmatrix, TdebugPolicy>"
176  "::step(): HOMO-LUMO gap has not increased in the "
177  "last five iterations. Desired accuracy too high for"
178  "current floating point precision?");
179  }
180 
181  ++niter;
182  currentStep.setTimeTotal(tStepTotal.toc());
183  }
184 
185  template<typename Treal, typename Tmatrix, typename TdebugPolicy>
187 
188  while (niter < info.getMaxSteps() - 1 && !info.converged() ) {
189  step();
190  }
191  if ( info.converged() ) {
192  // Only if converged because forcing correct occupation can have
193  // strange effects otherwise.
194  info.forceCorrectOccupation();
195  info.improveCorrectOccupation();
196  info.improveInfo();
197  info.improveHomoLumoF();
198  }
199  }
200 
201  template<typename Treal, typename Tmatrix, typename TdebugPolicy>
204  Treal const XmX2ENIsSmallValue = 0.207106781186547;
205  Treal const ONE = 1.0;
206 
207  Time t;
208  t.tic();
209  X2 = ONE * X * X;
210  currentStep.setTimeSquare(t.toc());
211 
212  currentStep.setNnzX(X.nnz());
213  currentStep.setNnzX2(X2.nnz());
214  currentStep.computeEigAccLoss();
215  currentStep.computeExactValues(X, X2);
216  currentStep.setTraceX(X.trace());
217  currentStep.setTraceX2(X2.trace());
218 
219  /* Now we are about to compute the Euclidean norm of (X -
220  X2). Previously this operation needed lots of memory. Now,
221  however, we hope that the memory usage is smaller because the
222  difference matrix is never explicitly computed. */
223 
224  typename Tmatrix::VectorType * eigVecPtr = new typename Tmatrix::VectorType;
225  Treal diffAcc;
226  Interval<Treal> XmX2EN;
227  t.tic();
228  if (info.ShouldComputeXmX2EuclNormAccurately(diffAcc)) {
229  if (normXmX2 == euclNorm)
230  XmX2EN = Tmatrix::euclDiffIfSmall(X, X2,
231  diffAcc,
232  XmX2ENIsSmallValue,
233  eigVecPtr);
234  else
235  XmX2EN = Tmatrix::diffIfSmall(X, X2,
236  normXmX2, diffAcc,
237  XmX2ENIsSmallValue);
238  }
239  else {
240  XmX2EN = Tmatrix::diffIfSmall(X, X2,
241  frobNorm, diffAcc,
242  XmX2ENIsSmallValue);
243  }
244  currentStep.setTimeXmX2Norm(t.toc());
245 
246  if (!eigVecPtr->is_empty())
247  currentStep.setEigVecPtr(eigVecPtr);
248  else
249  delete eigVecPtr;
250 
251  XmX2EN.increase(currentStep.getEigAccLoss());
252  Interval<Treal> zeroInt(0.0, XmX2EN.upp());
253  XmX2EN.intersect(zeroInt);
254  // FIXME:
255  // Consider improving lanczos so that if to good accuracy is requested
256  // the best possible accuracy is returned
257  currentStep.setXmX2EuclNorm(XmX2EN);
258  currentStep.checkIntervals("Purification::stepComputeInfo after setXmX2EuclNorm.");
259 
260  Treal tmpVal = template_blas_sqrt(1 + 4 * XmX2EN.upp());
261  currentStep.improveEigInterval
262  (Interval<Treal>((1 - tmpVal) / 2, (1 + tmpVal) / 2));
263 
264  currentStep.checkIntervals("Purification::stepComputeInfo B");
265 
266  info.improveInfo();
267  if (currentStep.getEigInterval().length() > 1.5)
268  throw Failure("Purification<Treal, Tmatrix, TdebugPolicy>"
269  "::step() : "
270  "Eigenvalues moved to far from [0, 1] interval.");
271  /* Now decide which polynomial to use for first step. */
272  currentStep.setPoly();
273  currentStep.checkIntervals("Purification::stepComputeInfo C");
274  }
275 
276 
277 
278 
279 } /* end namespace mat */
280 #undef __ID__
281 #endif
void setActualThresh(Treal const thr)
Definition: PuriStepInfo.h:85
Tmatrix & X
Definition: Purification.h:61
Treal getEigAccLoss() const
Definition: PuriStepInfo.h:140
Interval< Treal > getHomoF() const
Returns the best interval containing the homo eigenvalue (not transformed to [0, 1]) ...
Definition: PuriInfo.h:112
Treal upp() const
Definition: Interval.h:143
void computeExactValues(Tmatrix const &X, Tmatrix const &X2)
Definition: PuriStepInfo.h:152
void purify()
Definition: Purification.h:186
void setTimeThresh(float t)
Definition: PuriStepInfo.h:164
void improveEigInterval(Interval< Treal > const eInt)
Improve eigenvalue bounds and delta if possible.
Definition: PuriStepInfo.h:436
void tic()
Definition: matInclude.cc:77
static Interval intersect(Interval const &A, Interval const &B)
Definition: Interval.h:51
Interval< Treal > getLumoF() const
Returns the best interval containing the lumo eigenvalue (not transformed to [0, 1]) ...
Definition: PuriInfo.h:116
PuriInfo< Treal, VectorType, TdebugPolicy > & info
Definition: Purification.h:65
Interval< Treal > const & getLumo() const
Definition: PuriStepInfo.h:118
void improveHomoLumo(Interval< Treal > const homoInt, Interval< Treal > const lumoInt)
Definition: PuriStepInfo.h:310
void increase(Treal const value)
Increases interval with value in both directions.
Definition: Interval.h:131
bool empty() const
Definition: Interval.h:49
void setTimeXmX2Norm(float t)
Definition: PuriStepInfo.h:166
Definition: matInclude.h:158
float toc()
Definition: matInclude.cc:81
Definition: Interval.h:44
void setTimeTotal(float t)
Definition: PuriStepInfo.h:167
Definition: matInclude.h:135
void setChosenThresh(Treal const thr)
Definition: PuriStepInfo.h:83
void setTraceX(Treal const trX)
Definition: PuriStepInfo.h:91
void setTimeSquare(float t)
Definition: PuriStepInfo.h:165
Treal low() const
Definition: Interval.h:142
Tmatrix::VectorType VectorType
Definition: Purification.h:45
Definition: Purification.h:43
PuriStepInfo< Treal, Tvector, TdebugPolicy > & getNext()
Definition: PuriInfo.h:84
Interval< Treal > const & getHomo() const
Definition: PuriStepInfo.h:115
normType const normXmX2
Definition: Purification.h:64
void stepComputeInfo(PuriStepInfo< Treal, VectorType, TdebugPolicy > &currentStep)
Definition: Purification.h:203
void setTraceX2(Treal const trX2)
Definition: PuriStepInfo.h:92
void setNnzX(size_t const nzX)
Definition: PuriStepInfo.h:135
void setPoly()
Definition: PuriStepInfo.h:292
void setEigVecPtr(Tvector *eigVecPtr_)
Improves homo and lumo bounds if the new ones are better.
Definition: PuriStepInfo.h:107
Treal getOptimalThresh() const
Definition: PuriInfo.h:363
Definition: matInclude.h:135
Tmatrix X2
Definition: Purification.h:62
Definition: Failure.h:47
void step()
Definition: Purification.h:121
void checkIntervals(const char *descriptionString) const
Definition: PuriStepInfo.h:147
void setNnzX2(size_t const nzX2)
Definition: PuriStepInfo.h:137
int niter
Definition: Purification.h:66
Purification(Tmatrix &M, normType const normXmX2, PuriInfo< Treal, VectorType, TdebugPolicy > &info)
Constructor.
Definition: Purification.h:75
void computeEigAccLoss()
Computes a probable upper bound of the accuracy that is lost in the eigenvalues of X * X because of l...
Definition: PuriStepInfo.h:467
void setXmX2EuclNorm(Interval< Treal > const XmX2eucl)
Sets XmX2EuclNorm bounds.
Definition: PuriStepInfo.h:98
Interval< Treal > const & getEigInterval() const
Definition: PuriStepInfo.h:121
void setMemUsageBeforeTrunc()
Definition: PuriStepInfo.h:158
Treal template_blas_sqrt(Treal x)
normType const normTruncation
Definition: Purification.h:63
normType
Definition: matInclude.h:135
Interval< Treal > getEigFInterval() const
Definition: PuriInfo.h:81