ergo
PuriStepInfoDebug.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_PURISTEPINFODEBUG
37 #define MAT_PURISTEPINFODEBUG
38 #include "matInclude.h"
39 #include "Interval.h"
40 #define __ID__ "$Id$"
41 namespace mat {
42  template<typename Treal, typename TdebugPolicy>
43  class PuriStepInfoDebug : public TdebugPolicy {
44  public:
45 
46  inline void checkIntervals(Interval<Treal> const & eigInterval,
47  Interval<Treal> const & homo,
48  Interval<Treal> const & lumo,
49  Interval<Treal> const & XmX2EuclNorm,
50  const char* descriptionString) const {}
51  template<typename Tmatrix>
52  inline void computeExactValues(Tmatrix const & X, Tmatrix const & X2,
53  int const n, int const nocc) {}
54  protected:
56  private:
57  };
58 
59 
60 
61 /* Specialization for high debug level */
62  template<typename Treal>
64  public:
65  void checkIntervals(Interval<Treal> const & eigInterval,
66  Interval<Treal> const & homo,
67  Interval<Treal> const & lumo,
68  Interval<Treal> const & XmX2EuclNorm,
69  const char* descriptionString) const;
70  template<typename Tmatrix>
71  void computeExactValues(Tmatrix const & X, Tmatrix const & X2,
72  int const n, int const nocc);
73  protected:
75  :homoExact(), lumoExact(),
76  lmaxExact(), lminExact(),
77  XmX2EuclNormExact()
78  {}
84  private:
85  };
86 
87  template<typename Treal>
89  checkIntervals(Interval<Treal> const & eigInterval,
90  Interval<Treal> const & homo,
91  Interval<Treal> const & lumo,
92  Interval<Treal> const & XmX2EuclNorm,
93  const char* descriptionString) const {
94  /* Failure here probably means that checkIntervals was called before
95  * the exact values were calculated.
96  */
97  if(lminExact.empty())
98  std::cout << "PuriStepInfoDebug::checkIntervals failed, descriptionString = '"
99  << descriptionString << "'" << std::endl;
100  ASSERTALWAYS(!lminExact.empty());
101  Interval<Treal> zeroOneInt(0.0,1.0);
102  if (Interval<Treal>::intersect(eigInterval, lminExact).empty())
103  std::cout<<std::endl<<" eigInterval "<<
104  std::setprecision(25)<<eigInterval<<std::endl
105  <<" lminExact "<<std::setprecision(25)
106  <<lminExact<<std::endl;
107  ASSERTDEBUG(!Interval<Treal>::intersect(eigInterval, lminExact).empty());
108  if (Interval<Treal>::intersect(eigInterval, lmaxExact).empty())
109  std::cout<<std::endl<<" eigInterval "<<
110  std::setprecision(25)<<eigInterval<<std::endl
111  <<" lmaxExact "<<lmaxExact<<std::endl;
112  ASSERTDEBUG(!Interval<Treal>::intersect(eigInterval, lmaxExact).empty());
113  /* The homo and lumo intervals only provides information if
114  * the exact values lie in [0, 1] otherwise no check is done.
115  */
116  if (!Interval<Treal>::intersect(zeroOneInt, homoExact).empty()) {
117  if (Interval<Treal>::intersect(homo, homoExact).empty())
118  std::cout<<std::endl<<" homo "<<
119  std::setprecision(25)<<homo<<std::endl
120  <<" homoExact "<<homoExact<<std::endl;
121  if(Interval<Treal>::intersect(homo, homoExact).empty())
122  std::cout << "PuriStepInfoDebug::checkIntervals failed due to intersect(homo, homoExact) , descriptionString = '"
123  << descriptionString << "'" << std::endl;
124  ASSERTDEBUG(!Interval<Treal>::intersect(homo, homoExact).empty());
125  }
126  if (!Interval<Treal>::intersect(zeroOneInt, lumoExact).empty()) {
127  if (Interval<Treal>::intersect(lumo, lumoExact).empty())
128  std::cout<<std::endl<<" lumo "<<
129  std::setprecision(25)<<lumo<<std::endl
130  <<" lumoExact "<<lumoExact<<std::endl;
131  ASSERTDEBUG(!Interval<Treal>::intersect(lumo, lumoExact).empty());
132  }
133  if (Interval<Treal>::
134  intersect(XmX2EuclNorm, XmX2EuclNormExact).empty()) {
135  std::cout<<std::endl<<" XmX2EuclNorm "<<
136  std::setprecision(25)<<XmX2EuclNorm<<std::endl
137  <<" XmX2EuclNormExact "<<XmX2EuclNormExact<<std::endl;
138  std::cout << "PuriStepInfoDebug::checkIntervals failed, descriptionString = '"
139  << descriptionString << "'" << std::endl;
140  }
142  intersect(XmX2EuclNorm, XmX2EuclNormExact).empty());
143  }
144 
145  template<typename Treal>
146  template<typename Tmatrix>
148  computeExactValues(Tmatrix const & X, Tmatrix const & X2,
149  int const n, int const nocc) {
150  /* Set up work space */
151  std::vector<Treal> full(n*n);
152  // Treal* full = new Treal[n*n];
153  Treal* eigvals = new Treal[n];
154  int lwork = 3*n-1;
155  Treal* work = new Treal[lwork];
156  int info = 0;
157  Treal euclNormMatrix;
158  Treal precSyev;
159  /* Compute lumoExact, homoExact, lminExact, and lmaxExact. */
160  X.fullMatrix(full);
161  syev("N", "U", &n, &full[0], &n, eigvals, work, &lwork, &info);
162  ASSERTALWAYS(!info);
163  euclNormMatrix = template_blas_fabs(eigvals[0]) > template_blas_fabs(eigvals[n-1]) ?
164  template_blas_fabs(eigvals[0]) : template_blas_fabs(eigvals[n-1]);
165  precSyev = euclNormMatrix * getRelPrecision<Treal>();
166  lumoExact = Interval<Treal>(eigvals[n-nocc-1] - precSyev,
167  eigvals[n-nocc-1] + precSyev);
168  homoExact = Interval<Treal>(eigvals[n-nocc] - precSyev,
169  eigvals[n-nocc] + precSyev);
170  lminExact = Interval<Treal>(eigvals[0] - precSyev,
171  eigvals[0] + precSyev);
172  lmaxExact = Interval<Treal>(eigvals[n-1] - precSyev,
173  eigvals[n-1] + precSyev);
174  /* Compute largest magnitude eigenvalue of X-X^2 */
175  Tmatrix XmX2(X);
176  XmX2 += -1.0 * X2;
177  XmX2.fullMatrix(full);
178  syev("N", "U", &n, &full[0], &n, eigvals, work, &lwork, &info);
179  ASSERTALWAYS(!info);
180  euclNormMatrix = template_blas_fabs(eigvals[0]) > template_blas_fabs(eigvals[n-1]) ?
181  template_blas_fabs(eigvals[0]) : template_blas_fabs(eigvals[n-1]);
182  precSyev = euclNormMatrix * getRelPrecision<Treal>();
183 
184 #if 0
185  std::cout<<"EIGS XmX2: "<<std::endl;
186  for(int ind = 0; ind < n; ++ind)
187  std::cout<<std::setprecision(20)<<eigvals[ind]<<std::endl;
188  std::cout<<"end EIGS XmX2: "<<std::endl<<std::endl;
189 #endif
190 
191  Treal lowBound = euclNormMatrix - precSyev > 0 ?
192  euclNormMatrix - precSyev : 0;
193  XmX2EuclNormExact = Interval<Treal>(lowBound, euclNormMatrix + precSyev);
194  /* Free memory */
195  delete[] eigvals;
196  delete[] work;
197  }
198 
199 
200 } /* end namespace mat */
201 #undef __ID__
202 #endif
Interval< Treal > XmX2EuclNormExact
Definition: PuriStepInfoDebug.h:83
Definition: PuriStepInfoDebug.h:43
PuriStepInfoDebug()
Definition: PuriStepInfoDebug.h:74
Interval< Treal > lminExact
Definition: PuriStepInfoDebug.h:82
Interval< Treal > lmaxExact
Definition: PuriStepInfoDebug.h:81
#define ASSERTALWAYS(x)
Definition: DebugPolicies.h:83
Treal template_blas_fabs(Treal x)
Definition: Interval.h:44
Interval< Treal > homoExact
Definition: PuriStepInfoDebug.h:79
static void syev(const char *jobz, const char *uplo, const int *n, T *a, const int *lda, T *w, T *work, const int *lwork, int *info)
Definition: gblas.h:380
PuriStepInfoDebug()
Definition: PuriStepInfoDebug.h:55
#define ASSERTDEBUG(x)
Definition: DebugPolicies.h:85
void checkIntervals(Interval< Treal > const &eigInterval, Interval< Treal > const &homo, Interval< Treal > const &lumo, Interval< Treal > const &XmX2EuclNorm, const char *descriptionString) const
Definition: PuriStepInfoDebug.h:46
Definition: DebugPolicies.h:89
Copyright(c) Emanuel Rubensson 2006.
void computeExactValues(Tmatrix const &X, Tmatrix const &X2, int const n, int const nocc)
Definition: PuriStepInfoDebug.h:52
Interval class.
Interval< Treal > lumoExact
Definition: PuriStepInfoDebug.h:80