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