ergo
Perturbation.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_PERTURBATION
37 #define MAT_PERTURBATION
38 namespace per {
39  template<typename Treal, typename Tmatrix, typename Tvector>
40  class Perturbation {
41  public:
42  Perturbation(std::vector<Tmatrix *> const & F,
44  std::vector<Tmatrix *> & D,
46  mat::Interval<Treal> const & gap,
53  Treal const deltaMax,
54  Treal const errorTol,
55  mat::normType const norm,
56  Tvector & vect
57  );
58  void perturb() {
59  dryRun();
60  run();
61  }
62 
63  void checkIdempotencies(std::vector<Treal> & idemErrors);
64  template<typename TmatNoSymm>
65  void checkCommutators(std::vector<Treal> & commErrors,
66  TmatNoSymm const & dummyMat);
67  void checkMaxSubspaceError(Treal & subsError);
68 
69  protected:
70  /* This is input from the beginning */
71  std::vector<Tmatrix *> const & F;
72  std::vector<Tmatrix *> & X;
75  Treal deltaMax;
76  Treal errorTol;
78  Tvector & vect;
79 
80  /* These variables are set in the dry run. */
81  int nIter;
82  std::vector<Treal> threshVal;
83  std::vector<Treal> sigma;
84 
95  void dryRun();
96  void run();
97  private:
98 
99  };
100 
101  template<typename Treal, typename Tmatrix, typename Tvector>
103  Perturbation(std::vector<Tmatrix *> const & F_in,
104  std::vector<Tmatrix *> & X_in,
105  mat::Interval<Treal> const & gap_in,
106  mat::Interval<Treal> const & allEigs_in,
107  Treal const deltaMax_in,
108  Treal const errorTol_in,
109  mat::normType const norm_in,
110  Tvector & vect_in)
111  : F(F_in), X(X_in), gap(gap_in), allEigs(allEigs_in),
112  deltaMax(deltaMax_in), errorTol(errorTol_in), norm(norm_in),
113  vect(vect_in) {
114  if (!X.empty())
115  throw "Perturbation constructor: D vector is expected to be empty (size==0)";
116  for (unsigned int iMat = 0; iMat < F.size(); ++iMat)
117  X.push_back(new Tmatrix(*F[iMat]));
118 
119  Treal lmin = allEigs.low();
120  Treal lmax = allEigs.upp();
121 
122  /***** Initial linear transformation of matrix sequence. */
123  typename std::vector<Tmatrix *>::iterator matIt = X.begin();
124  /* Scale to [0, 1] interval and negate */
125  (*matIt)->add_identity(-lmax);
126  *(*matIt) *= ((Treal)1.0 / (lmin - lmax));
127  matIt++;
128  /* ...and derivatives: */
129  for ( ; matIt != X.end(); matIt++ )
130  *(*matIt) *= ((Treal)-1.0 / (lmin - lmax));
131  /* Compute transformed gap */
132  gap = (gap - lmax) / (lmin - lmax);
133  }
134 
135  template<typename Treal, typename Tmatrix, typename Tvector>
137  Treal errorTolPerIter;
138  int nIterGuess = 0;
139  nIter = 1;
140  Treal lumo;
141  Treal homo;
142  Treal m;
143  Treal g;
144  while (nIterGuess < nIter) {
145  nIterGuess++;
146  errorTolPerIter = 0.5 * errorTol /nIterGuess;
147  nIter = 0;
148  mat::Interval<Treal> gapTmp(gap);
149  sigma.resize(0);
150  threshVal.resize(0);
151  while (gapTmp.low() > 0.5 * errorTol || gapTmp.upp() < 0.5 * errorTol) {
152  lumo = gapTmp.low();
153  homo = gapTmp.upp();
154  m = gapTmp.midPoint();
155  g = gapTmp.length();
156  if (m > 0.5) {
157  lumo = lumo*lumo;
158  homo = homo*homo;
159  sigma.push_back(-1);
160  }
161  else {
162  lumo = 2*lumo - lumo*lumo;
163  homo = 2*homo - homo*homo;
164  sigma.push_back(1);
165  }
166  /* Compute threshold value necessary to converge. */
167  Treal forceConvThresh = template_blas_fabs(1-2*m) * g / 10;
168  /* We divide by 10 > 2 so that this loop converges at some point. */
169  /* Compute threshold value necessary to maintain accuracy in subspace.*/
170  Treal subspaceThresh = errorTolPerIter * (homo-lumo) / (1+errorTolPerIter);
171  /* Choose the most restrictive threshold of the two. */
172  threshVal.push_back(forceConvThresh < subspaceThresh ?
173  forceConvThresh : subspaceThresh);
174  homo -= threshVal.back();
175  lumo += threshVal.back();
176  gapTmp = mat::Interval<Treal>(lumo, homo);
177  if (gapTmp.empty())
178  throw "Perturbation<Treal, Tmatrix, Tvector>::dryRun() : Perturbation iterations will fail to converge; Gap is too small or desired accuracy too high.";
179  nIter++;
180  } /* end 2nd while */
181  } /* end 1st while */
182  /* Now, we have nIter, threshVal, and sigma. */
183  }
184 
185  template<typename Treal, typename Tmatrix, typename Tvector>
187  Treal const ONE = 1.0;
188  mat::SizesAndBlocks rowsCopy;
189  X.front()->getRows(rowsCopy);
190  mat::SizesAndBlocks colsCopy;
191  X.front()->getCols(colsCopy);
192  Tmatrix tmpMat;
193  // tmpMat.resetSizesAndBlocks(rowsCopy, colsCopy);
194  int nMatrices;
195  Treal threshValPerOrder;
196  Treal chosenThresh;
197  for (int iter = 0; iter < nIter; iter++) {
198  std::cout<<"\n\nInside outer loop iter = "<<iter
199  <<" nIter = "<<nIter
200  <<" sigma = "<<sigma[iter]<<std::endl;
201  /* Number of matrices increases by 1 in each iteration: */
202  X.push_back(new Tmatrix);
203  nMatrices = X.size();
204  X[nMatrices-1]->resetSizesAndBlocks(rowsCopy, colsCopy);
205  /* Compute threshold value for each order. */
206  threshValPerOrder = threshVal[iter] / nMatrices;
207  /* Loop through all nonzero orders. */
208  std::cout<<"Entering inner loop nMatrices = "<<nMatrices<<std::endl;
209  for (int j = nMatrices-1 ; j >= 0 ; --j) {
210  std::cout<<"Inside inner loop j = "<<j<<std::endl;
211  std::cout<<"X[j]->eucl() (before compute) = "<<X[j]->eucl(vect,1e-7)<<std::endl;
212  std::cout<<"X[j]->frob() (before compute) = "<<X[j]->frob()<<std::endl;
213  tmpMat = Treal(Treal(1.0)+sigma[iter]) * (*X[j]);
214  std::cout<<"tmpMat.eucl() (before for) = "<<tmpMat.eucl(vect,1e-7)<<std::endl;
215  std::cout<<"tmpMat.frob() (before for) = "<<tmpMat.frob()<<std::endl;
216  for (int k = 0; k <= j; k++) {
217  /* X[j] = X[j] - sigma * X[k] * X[j-k] */
218  Tmatrix::ssmmUpperTriangleOnly(-sigma[iter], *X[k], *X[j-k],
219  ONE, tmpMat);
220  } /* End 3rd for */
221  std::cout<<"tmpMat.eucl() (after for) = "<<tmpMat.eucl(vect,1e-7)<<std::endl;
222  *X[j] = tmpMat;
223 
224  /* Truncate tmpMat, remove if it becomes zero. */
225  chosenThresh = threshValPerOrder / pow(deltaMax, Treal(j));
226  std::cout<<"X[j]->eucl() (before thresh) = "<<X[j]->eucl(vect,1e-7)<<std::endl;
227  std::cout<<"Chosen thresh: "<<chosenThresh<<std::endl;
228  Treal actualThresh = X[j]->thresh(chosenThresh, vect, norm);
229  std::cout<<"X[j]->eucl() (after thresh) = "<<X[j]->eucl(vect,1e-7)<<std::endl;
230 #if 1
231  /* If the current matrix is zero AND
232  * it is the last matrix
233  */
234  if (*X[j] == 0 && int(X.size()) == j+1) {
235  std::cout<<"DELETION: j = "<<j<<" X.size() = "<<X.size()
236  <<" X[j] = "<<X[j]<< " X[j]->frob() = "<<X[j]->frob()
237  <<std::endl;
238  delete X[j];
239  X.pop_back();
240  }
241  else
242  std::cout<<"NO DELETION: j = "<<j<<" X.size() = "<<X.size()
243  <<" X[j] = "<<X[j]<< " X[j]->frob() = "<<X[j]->frob()
244  <<std::endl;
245 #endif
246 
247  } /* End 2nd for (Loop through orders) */
248  } /* End 1st for (Loop through iterations) */
249  } /* End run() */
250 
251 
252  template<typename Treal, typename Tmatrix, typename Tvector>
254  checkIdempotencies(std::vector<Treal> & idemErrors) {
255  Tmatrix tmpMat;
256  Treal const ONE = 1.0;
257  unsigned int j;
258  for (unsigned int m = 0; m < X.size(); ++m) {
259  tmpMat = (-ONE) * (*X[m]);
260  for (unsigned int i = 0; i <= m; ++i) {
261  j = m - i;
262  /* TMP = TMP + X[i] * X[j] */
263  Tmatrix::ssmmUpperTriangleOnly(ONE, *X[i], *X[j], ONE, tmpMat);
264  }
265  /* TMP is symmetric! */
266  idemErrors.push_back(tmpMat.eucl(vect,1e-10));
267  }
268  }
269 
270  template<typename Treal, typename Tmatrix, typename Tvector>
271  template<typename TmatNoSymm>
273  checkCommutators(std::vector<Treal> & commErrors,
274  TmatNoSymm const & dummyMat) {
275  mat::SizesAndBlocks rowsCopy;
276  X.front()->getRows(rowsCopy);
277  mat::SizesAndBlocks colsCopy;
278  X.front()->getCols(colsCopy);
279  TmatNoSymm tmpMat;
280  tmpMat.resetSizesAndBlocks(rowsCopy, colsCopy);
281  Treal const ONE = 1.0;
282  unsigned int j;
283  for (unsigned int m = 0; m < X.size(); ++m) {
284  tmpMat = 0;
285  std::cout<<"New loop\n";
286  for (unsigned int i = 0; i <= m && i < F.size(); ++i) {
287  j = m - i;
288  std::cout<<i<<", "<<j<<std::endl;
289  /* TMP = TMP + F[i] * X[j] - X[j] * F[i] */
290  tmpMat += ONE * (*F[i]) * (*X[j]);
291  tmpMat += -ONE * (*X[j]) * (*F[i]);
292  }
293  /* TMP is not symmetric! */
294  commErrors.push_back(tmpMat.frob());
295  }
296  }
297 
298 
299  template<typename Treal, typename Tmatrix, typename Tvector>
301  checkMaxSubspaceError(Treal & subsError) {
302  Treal const ONE = 1.0;
303  Tmatrix XdeltaMax(*F.front());
304  for (unsigned int ind = 1; ind < F.size(); ++ind)
305  XdeltaMax += pow(deltaMax, Treal(ind)) * (*F[ind]);
306  /***** Initial linear transformation of matrix. */
307  Treal lmin = allEigs.low();
308  Treal lmax = allEigs.upp();
309  /* Scale to [0, 1] interval and negate */
310  XdeltaMax.add_identity(-lmax);
311  XdeltaMax *= ((Treal)1.0 / (lmin - lmax));
312 
313  Tmatrix X2;
314  for (int iter = 0; iter < nIter; iter++) {
315  X2 = ONE * XdeltaMax * XdeltaMax;
316  if (sigma[iter] == Treal(1.0)) {
317  XdeltaMax *= 2.0;
318  XdeltaMax -= X2;
319  }
320  else {
321  XdeltaMax = X2;
322  }
323  } /* End of for (Loop through iterations) */
324 
325  Tmatrix DdeltaMax(*X.front());
326  for (unsigned int ind = 1; ind < X.size(); ++ind)
327  DdeltaMax += pow(deltaMax, Treal(ind)) * (*X[ind]);
328  subsError = Tmatrix::eucl_diff(XdeltaMax,DdeltaMax,
329  vect, errorTol *1e-2);
330  }
331 
332 
333 
334 } /* end namespace mat */
335 #endif
336 
Tvector & vect
Definition: Perturbation.h:78
void dryRun()
Dry run to obtain some needed numbers.
Definition: Perturbation.h:136
Treal upp() const
Definition: Interval.h:143
Treal length() const
Returns the length of the interval.
Definition: Interval.h:107
std::vector< Tmatrix * > const & F
Definition: Perturbation.h:71
Perturbation(std::vector< Tmatrix * > const &F, std::vector< Tmatrix * > &D, mat::Interval< Treal > const &gap, mat::Interval< Treal > const &allEigs, Treal const deltaMax, Treal const errorTol, mat::normType const norm, Tvector &vect)
Definition: Perturbation.h:103
mat::Interval< Treal > const & allEigs
Definition: Perturbation.h:74
void run()
Definition: Perturbation.h:186
bool empty() const
Definition: Interval.h:49
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:37
Treal midPoint() const
Definition: Interval.h:113
Treal template_blas_fabs(Treal x)
Definition: Interval.h:44
void checkCommutators(std::vector< Treal > &commErrors, TmatNoSymm const &dummyMat)
Definition: Perturbation.h:273
Treal errorTol
Definition: Perturbation.h:76
void perturb()
Definition: Perturbation.h:58
std::vector< Tmatrix * > & X
Definition: Perturbation.h:72
void checkMaxSubspaceError(Treal &subsError)
Definition: Perturbation.h:301
mat::normType const norm
Definition: Perturbation.h:77
Treal low() const
Definition: Interval.h:142
Definition: Perturbation.h:40
std::vector< Treal > sigma
Definition: Perturbation.h:83
void checkIdempotencies(std::vector< Treal > &idemErrors)
Definition: Perturbation.h:254
Treal deltaMax
Definition: Perturbation.h:75
int nIter
Definition: Perturbation.h:81
std::vector< Treal > threshVal
Definition: Perturbation.h:82
mat::Interval< Treal > gap
Definition: Perturbation.h:73
normType
Definition: matInclude.h:135