ergo
purification_old.h
Go to the documentation of this file.
1 /* Ergo, version 3.4, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2014 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 
34 #ifndef MAT_PURIFICATION_OLD
35 #define MAT_PURIFICATION_OLD
36 #include <iomanip>
37 #include <math.h>
38 namespace mat{
39 
40 
41  template<class Treal, class MAT>
42  void tc2_extra(MAT& D, const int nshells,
43  const int nsteps, const Treal frob_trunc = 0) {
44  MAT tmp(D);
45  for (int step = 0; step < nsteps; step++) {
46  Treal tracediff = tmp.trace() - nshells;
47  if (tracediff > 0) {
48  tmp = (Treal)1.0 * D * D;
49  D = tmp;
50  D = (Treal)-1.0 * tmp * tmp + (Treal)2.0 * D;
51  }
52  else {
53  tmp = D;
54  tmp = (Treal)-1.0 * D * D + (Treal)2.0 * tmp;
55  D = (Treal)1.0 * tmp * tmp;
56  }
57  D.frob_thresh(frob_trunc);
58  }
59  }
60 
61  template<class Treal, class MAT>
62  void tc2_extra_auto(MAT& D, const int nshells,
63  int& nsteps, const Treal frob_trunc,
64  const int maxiter = 50) {
65  MAT tmp(D);
66  Treal froberror;
67  Treal tracediff = tmp.trace() - nshells;
68  nsteps = 0;
69  do {
70  if (tracediff > 0) {
71  tmp = (Treal)1.0 * D * D;
72  D = tmp;
73  D = (Treal)-1.0 * tmp * tmp + (Treal)2.0 * D;
74  }
75  else {
76  tmp = D;
77  tmp = (Treal)-1.0 * D * D + (Treal)2.0 * tmp;
78  D = (Treal)1.0 * tmp * tmp;
79  }
80  froberror = MAT::frob_diff(D, tmp);
81  D.frob_thresh(frob_trunc);
82  tracediff = D.trace() - nshells;
83 #if 0
84  std::cout<<"Iteration:"<<nsteps<<" Froberror: "
85  <<std::setprecision(8)<<froberror<<
86  " Tracediff: "<<tracediff<<std::endl;
87 #endif
88  nsteps++;
89  if (nsteps > maxiter)
90  throw AcceptableMaxIter("Extra Auto Purification reached maxiter"
91  " without convergence", maxiter);
92  } while (froberror > frob_trunc);
93  }
94 
95 
96 
97  template<class Treal, class MAT>
98  void tc2(MAT& D, int& iterations,
99  const MAT& H, const int nshells,
100  const Treal trace_error_limit = 1e-2,
101  const Treal frob_trunc = 0,
102  int* polys= NULL,
103  const int maxiter = 100) {
104  MAT tmp(H);
105  D = H;
106  iterations = 0;
107  Treal tracediff = tmp.trace() - nshells;
108  Treal tracediff_old = 2.0 * trace_error_limit;
109 
110  while (template_blas_fabs(tracediff) > trace_error_limit &&
111  template_blas_fabs(tracediff_old) > trace_error_limit) {
112  // std::cout<<"Iteration:"<<iterations
113  // <<" Tracediff ="<<tracediff<<std::endl;
114  if (tracediff > 0) {
115  D = (Treal)1.0 * tmp * tmp;
116  if (polys)
117  polys[iterations] = 0;
118  }
119  else {
120  D = (Treal)-1.0 * tmp * tmp + (Treal)2.0 * D;
121  if (polys)
122  polys[iterations] = 1;
123  }
124  D.frob_thresh(frob_trunc);
125  tmp = D;
126  tracediff_old = tracediff;
127  tracediff = D.trace() - nshells;
128  iterations++;
129  if (iterations > maxiter)
130  throw AcceptableMaxIter("Purification reached maxiter without convergence", maxiter);
131  }
132  }
133 
134 #if 1
135  template<class Treal, class MAT>
136  void tc2_auto(MAT& D, int& iterations,
137  const MAT& H, const int nocc,
138  const Treal frob_trunc = 0,
139  int* polys= NULL,
140  const int maxiter = 100) {
141  assert(frob_trunc >= 0);
142  assert(nocc >= 0);
143  assert(maxiter >= 0);
144  MAT tmp(H);
145  D = H;
146  iterations = 0;
147  Treal tracediff = 2;
148  Treal newtracediff = tmp.trace() - nocc;
149  Treal froberror = 1;
150  Treal tracechange = 0; /* Initial value has no relevance */
151  /* Need check for too loose truncation? */
152  while (2 * froberror > (3 - template_blas_sqrt((Treal)5)) / 2 - frob_trunc ||
153  template_blas_fabs(tracediff) > 1 ||
154  template_blas_fabs(tracechange) > (1 - 2 * froberror) * (1 - template_blas_fabs(tracediff))) {
155  // std::cout<<"Iteration:"<<iterations
156  // <<" Tracediff ="<<tracediff<<std::endl;
157  tracediff = newtracediff;
158  if (tracediff > 0) {
159  D = (Treal)1.0 * tmp * tmp;
160  if (polys)
161  polys[iterations] = 0;
162  }
163  else {
164  D = (Treal)-1.0 * tmp * tmp + (Treal)2.0 * D;
165  if (polys)
166  polys[iterations] = 1;
167  }
168  D.frob_thresh(frob_trunc);
169  newtracediff = D.trace() - nocc;
170  tracechange = newtracediff - tracediff;
171  froberror = MAT::frob_diff(D, tmp);
172  tmp = D;
173 #if 0
174  std::cout<<"Iteration:"<<iterations<<" epsilon: "
175  <<std::setprecision(8)<<std::setw(12)<<froberror<<
176  " delta: "<<std::setw(12)<<tracediff<<
177  " tracechange: "<<std::setw(12)<<tracechange<<std::endl;
178 #endif
179  iterations++;
180  if (iterations > maxiter)
181  throw AcceptableMaxIter("tc2_auto::Purification reached maxiter"
182  " without convergence", maxiter);
183  }
184 
185  /* Take one step to make sure the eigenvalues stays in */
186  /* { [ 0 , 2 * epsilon [ , ] 1 - 2 * epsilon , 1] } */
187  if (tracediff > 0) {
188  D = (Treal)-1.0 * tmp * tmp + (Treal)2.0 * D;
189  if (polys)
190  polys[iterations] = 1;
191  }
192  else {
193  D = (Treal)1.0 * tmp * tmp;
194  if (polys)
195  polys[iterations] = 0;
196  }
197  iterations++;
198 
199  /* Use second order convergence polynomials to converge completely */
200  D.frob_thresh(frob_trunc);
201  tracediff = D.trace() - nocc;
202  do {
203  if (tracediff > 0) {
204  tmp = (Treal)1.0 * D * D;
205  D = tmp;
206  D = (Treal)-1.0 * tmp * tmp + (Treal)2.0 * D;
207  if (polys) {
208  polys[iterations] = 0;
209  polys[iterations + 1] = 1;
210  }
211  }
212  else {
213  tmp = D;
214  tmp = (Treal)-1.0 * D * D + (Treal)2.0 * tmp;
215  D = (Treal)1.0 * tmp * tmp;
216  if (polys) {
217  polys[iterations] = 1;
218  polys[iterations + 1] = 0;
219  }
220  }
221  froberror = MAT::frob_diff(D, tmp);
222  D.frob_thresh(frob_trunc);
223  tracediff = D.trace() - nocc;
224 #if 0
225  std::cout<<"Iteration:"<<iterations<<" Froberror: "
226  <<std::setprecision(8)<<froberror<<
227  " Tracediff: "<<tracediff<<std::endl;
228 #endif
229  iterations += 2;
230  if (iterations > maxiter)
231  throw AcceptableMaxIter("tc2_auto::Purification reached maxiter"
232  " in extra part without convergence", maxiter);
233  } while (froberror > frob_trunc);
234 
235 
236  }
237 
238 #endif
239 
240 
241 
242 #if 0
243  iterations = 0;
244  /* D = (F - lmax * I) / (lmin - lmax) */
245  D = F;
246  D.add_identity(-lmax);
247  D *= (1.0 / (lmin - lmax));
248  /*****/ clock_t start;
249  /*****/ float syrktime = 0;
250  /*****/ float threshtime = 0;
251  MAT tmp;
252  Treal tracediff;
253  for (int steps = 0; steps < nextra_steps; steps++) {
254  tmp = D;
255  tracediff = tmp.trace() - nshells;
256  while (template_blas_fabs(tracediff) > trace_error_limit) {
257  iterations++;
258  /*****/ start = clock();
259  if (tracediff > 0)
260  MAT::syrk('U', false, 1.0, tmp, 0.0, D);
261  else
262  MAT::syrk('U', false, -1.0, tmp, 2.0, D);
263  /*****/ syrktime += ((float)(clock()-start))/(CLOCKS_PER_SEC);
264  D.sym_to_nosym();
265  /*****/ start = clock();
266  D.frob_thresh(frob_trunc);
267  /*****/ threshtime += ((float)(clock()-start))/(CLOCKS_PER_SEC);
268  tmp = D;
269  tracediff = tmp.trace() - nshells;
270  } /* end while */
271 
272  /* extra steps */
273  iterations += 2;
274  if (tracediff > 0) {
275  /*****/ start = clock();
276  MAT::syrk('U', false, 1.0, tmp, 0.0, D);
277  /*****/ syrktime += ((float)(clock()-start))/(CLOCKS_PER_SEC);
278  D.sym_to_nosym();
279  tmp = D;
280  /*****/ start = clock();
281  MAT::syrk('U', false, -1.0, tmp, 2.0, D);
282  /*****/ syrktime += ((float)(clock()-start))/(CLOCKS_PER_SEC);
283  }
284  else {
285  /*****/ start = clock();
286  MAT::syrk('U', false, -1.0, tmp, 2.0, D);
287  /*****/ syrktime += ((float)(clock()-start))/(CLOCKS_PER_SEC);
288  D.sym_to_nosym();
289  tmp = D;
290  /*****/ start = clock();
291  MAT::syrk('U', false, 1.0, tmp, 0.0, D);
292  /*****/ syrktime += ((float)(clock()-start))/(CLOCKS_PER_SEC);
293  }
294  D.sym_to_nosym();
295  /*****/ start = clock();
296  D.frob_thresh(frob_trunc);
297  /*****/ threshtime += ((float)(clock()-start))/(CLOCKS_PER_SEC);
298  } /* end for */
300  std::cout<<std::endl<<" total syrktime in tcp ="<<std::setw(5)
301  <<syrktime<<std::endl;
302  std::cout<<" total threshtime in tcp ="<<std::setw(5)
303  <<threshtime<<std::endl;
304 
305  }
306 #endif
307 } /* end namespace mat */
308 
309 #endif
Definition: Failure.h:71
void tc2_extra(MAT &D, const int nshells, const int nsteps, const Treal frob_trunc=0)
Definition: purification_old.h:42
static void syrk(const char *uplo, const char *trans, const int *n, const int *k, const T *alpha, const T *A, const int *lda, const T *beta, T *C, const int *ldc)
Definition: mat_gblas.h:332
Definition: allocate.cc:30
void tc2(MAT &D, int &iterations, const MAT &H, const int nshells, const Treal trace_error_limit=1e-2, const Treal frob_trunc=0, int *polys=NULL, const int maxiter=100)
Definition: purification_old.h:98
Treal template_blas_fabs(Treal x)
void tc2_extra_auto(MAT &D, const int nshells, int &nsteps, const Treal frob_trunc, const int maxiter=50)
Definition: purification_old.h:62
void tc2_auto(MAT &D, int &iterations, const MAT &H, const int nocc, const Treal frob_trunc=0, int *polys=NULL, const int maxiter=100)
Definition: purification_old.h:136
Treal template_blas_sqrt(Treal x)