ergo
Interval.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_INTERVAL
37 #define MAT_INTERVAL
38 #include <math.h>
39 #include <iomanip>
40 namespace mat {
41 
42  /* FIXME represent interval as midpoint and length to get better precision */
43  template<typename Treal>
44  class Interval {
45  public:
46  explicit Interval(Treal low = 1, Treal upp = -1)
48  }
49  inline bool empty() const {return lowerBound > upperBound;}
50 
51  static Interval intersect(Interval const & A, Interval const & B) {
52  if (A.empty())
53  return A;
54  else if (B.empty())
55  return B;
56  else
57  return Interval(A.lowerBound > B.lowerBound ?
58  A.lowerBound : B.lowerBound,
59  A.upperBound < B.upperBound ?
60  A.upperBound : B.upperBound);
61  }
62  inline void intersect(Interval const & other) {
63  if (this->empty())
64  return;
65  else if (other.empty()) {
66  *this = other;
67  return;
68  }
69  else {
70  this->lowerBound = this->lowerBound > other.lowerBound ?
71  this->lowerBound : other.lowerBound;
72  this->upperBound = this->upperBound < other.upperBound ?
73  this->upperBound : other.upperBound;
74  return;
75  }
76  }
77 
78  inline void intersect_always_non_empty(Interval const & other) {
79  if (this->empty()) {
80  *this = other;
81  return;
82  }
83  if (other.empty())
84  return;
85 
86  Interval intersection = intersect(*this, other);
87  if ( !intersection.empty() ) {
88  *this = intersection;
89  return;
90  }
91  // Ok, the intersection is actually empty
92  Treal tmp_val;
93  if ( this->lowerBound > other.upperBound )
94  tmp_val = (this->lowerBound + other.upperBound) / 2;
95  else if ( this->upperBound < other.lowerBound )
96  tmp_val = (this->upperBound + other.lowerBound) / 2;
97  else
98  assert(0); // This should never happen!
99  this->lowerBound = tmp_val;
100  this->upperBound = tmp_val;
101  return;
102  }
103 
107  inline Treal length() const {
108  if (empty())
109  return 0.0;
110  else
111  return upperBound - lowerBound;
112  }
113  inline Treal midPoint() const {
114  assert(!empty());
115  return (upperBound + lowerBound) / 2;
116  }
117  inline bool cover(Treal const value) const {
118  if (empty())
119  return false;
120  else
121  return (value <= upperBound && value >= lowerBound);
122  }
123  inline bool overlap(Interval const & other) const {
124  Interval tmp = intersect(*this, other);
125  return !tmp.empty();
126  }
127 
131  inline void increase(Treal const value) {
132  if (empty())
133  throw Failure("Interval<Treal>::increase(Treal const) : "
134  "Attempt to increase empty interval.");
135  lowerBound -= value;
136  upperBound += value;
137  }
138  inline void decrease(Treal const value) {
139  lowerBound += value;
140  upperBound -= value;
141  }
142  inline Treal low() const {return lowerBound;}
143  inline Treal upp() const {return upperBound;}
144 
145 #if 0
146  inline Interval<Treal> operator*(Interval<Treal> const & other) const {
147  assert(lowerBound >= 0);
148  assert(other.lowerBound >= 0);
149  return Interval<Treal>(lowerBound * other.lowerBound,
150  upperBound * other.upperBound);
151  }
152 #endif
153 
154  inline Interval<Treal> operator*(Treal const & value) const {
155  if (value < 0)
156  return Interval<Treal>(upperBound * value, lowerBound * value);
157  else
158  return Interval<Treal>(lowerBound * value, upperBound * value);
159  }
160 
161  /* Minus operator well defined? */
162  inline Interval<Treal> operator-(Interval<Treal> const & other) const {
163  return *this + (-1.0 * other);
164  // return Interval<Treal>(lowerBound - other.lowerBound,
165  // upperBound - other.upperBound);
166  }
167 
168  inline Interval<Treal> operator+(Interval<Treal> const & other) const {
169  return Interval<Treal>(lowerBound + other.lowerBound,
170  upperBound + other.upperBound);
171  }
172  inline Interval<Treal> operator/(Treal const & value) const {
173  if (value < 0)
174  return Interval<Treal>(upperBound / value, lowerBound / value);
175  else
176  return Interval<Treal>(lowerBound / value, upperBound / value);
177  }
178  inline Interval<Treal> operator-(Treal const & value) const {
179  return Interval<Treal>(lowerBound - value, upperBound - value);
180  }
181  inline Interval<Treal> operator+(Treal const & value) const {
182  return Interval<Treal>(lowerBound + value, upperBound + value);
183  }
184 
185 
186 
187  void puriStep(int poly);
188  void invPuriStep(int poly);
189  // The two functions above really not needed now, just special
190  // case of functions below with alpha == 1
191  void puriStep(int poly, Treal alpha);
192  void invPuriStep(int poly, Treal alpha);
193  protected:
194  Treal lowerBound;
195  Treal upperBound;
196  private:
197  };
198 
199 #if 0
200  /* Minus operator is not well defined for intervals */
201  template<typename Treal>
202  inline Interval<Treal> operator-(Treal const value,
203  Interval<Treal> const & other) {
204  return Interval<Treal>(value - other.lowerBound,
205  value - other.upperBound);
206  }
207  template<typename Treal>
208  inline Interval<Treal> operator+(Treal const value,
209  Interval<Treal> const & other) {
210  return Interval<Treal>(value + other.lowerBound,
211  value + other.upperBound);
212  }
213 #endif
214 
215 #if 1
216  template<typename Treal>
218  if (other.empty())
219  return other;
220  else {
221  assert(other.low() >= 0);
222  return Interval<Treal>(template_blas_sqrt(other.low()), template_blas_sqrt(other.upp()));
223  }
224  }
225 #endif
226 
227  template<typename Treal>
228  void Interval<Treal>::puriStep(int poly) {
229  if (upperBound > 2.0 || lowerBound < -1.0)
230  throw Failure("Interval<Treal>::puriStep(int) : It is assumed here "
231  "that the interval I is within [-1.0, 2.0]");
232  Interval<Treal> oneInt(-1.0,1.0);
233  Interval<Treal> zeroInt(0.0,2.0);
234  /* Elias note 2010-09-12:
235  Sometimes the polynomial 2*x-x*x makes a non-empty interval in [0,1]
236  become empty because of rounding errors. For some reason this problem
237  showed up when using Intel compiler version 11.1 but not when using
238  version 10.1. Many test cases failed on Kalkyl when using
239  the compiler "icpc (ICC) 11.1 20100414".
240  Anyway, we know that the function f(x) = 2*x-x*x is growing
241  in the interval [0,1] so we know a non-empty interval inside [0,1] should
242  always stay non-empty. To avoid assertion failures in purification,
243  the code below now forces the interval to be non-empty if it became
244  empty due to rounding errors. */
245  bool nonEmptyIntervalInZeroToOne = false;
246  if(upperBound > lowerBound && lowerBound >= 0 && upperBound <= 1)
247  nonEmptyIntervalInZeroToOne = true;
248  if (poly) {
249  /* 2*x - x*x */
250  *this = Interval<Treal>::intersect(*this,oneInt);
251  // assert(upperBound <= 1.0);
252  upperBound = 2 * upperBound - upperBound * upperBound;
253  lowerBound = 2 * lowerBound - lowerBound * lowerBound;
254  if(nonEmptyIntervalInZeroToOne && upperBound < lowerBound) {
255  // Force interval to be non-empty.
256  Treal midPoint = (lowerBound + upperBound) / 2;
257  upperBound = lowerBound = midPoint;
258  }
259  }
260  else { /* x*x */
261  *this = Interval<Treal>::intersect(*this,zeroInt);
262  // assert(lowerBound >= 0.0);
263  upperBound = upperBound * upperBound;
264  lowerBound = lowerBound * lowerBound;
265  }
266  }
267 
268  template<typename Treal>
270  if (upperBound > 2.0 || lowerBound < -1.0)
271  throw Failure("Interval<Treal>::puriStep(int) : It is assumed here "
272  "that the interval I is within [-1.0, 1.0]");
273  Interval<Treal> oneInt(-1.0,1.0);
274  Interval<Treal> zeroInt(0.0,2.0);
275  if (poly) {
276  /* 1 - sqrt(1 - x) */
277  *this = Interval<Treal>::intersect(*this,oneInt);
278  // assert(upperBound <= 1.0);
279  upperBound = 1.0 - template_blas_sqrt(1.0 - upperBound);
280  lowerBound = 1.0 - template_blas_sqrt(1.0 - lowerBound);
281  }
282  else { /* sqrt(x) */
283  *this = Interval<Treal>::intersect(*this,zeroInt);
284  // assert(lowerBound >= 0.0);
285  upperBound = template_blas_sqrt(upperBound);
286  lowerBound = template_blas_sqrt(lowerBound);
287  }
288  }
289 
290 #if 1
291  template<typename Treal>
292  void Interval<Treal>::puriStep(int poly, Treal alpha) {
293  if (upperBound > 2.0 || lowerBound < -1.0)
294  throw Failure("Interval<Treal>::puriStep(int, real) : It is assumed here "
295  "that the interval I is within [-1.0, 2.0]");
296  Interval<Treal> oneInt(-1.0,1.0);
297  Interval<Treal> zeroInt(0.0,2.0);
298  /* Elias note 2010-09-12:
299  Sometimes the polynomial 2*x-x*x makes a non-empty interval in [0,1]
300  become empty because of rounding errors. For some reason this problem
301  showed up when using Intel compiler version 11.1 but not when using
302  version 10.1. Many test cases failed on Kalkyl when using
303  the compiler "icpc (ICC) 11.1 20100414".
304  Anyway, we know that the function f(x) = 2*x-x*x is growing
305  in the interval [0,1] so we know a non-empty interval inside [0,1] should
306  always stay non-empty. To avoid assertion failures in purification,
307  the code below now forces the interval to be non-empty if it became
308  empty due to rounding errors. */
309  bool nonEmptyIntervalInZeroToOne = false;
310  if(upperBound > lowerBound && lowerBound >= 0 && upperBound <= 1)
311  nonEmptyIntervalInZeroToOne = true;
312  if (poly) {
313  /* 2*alpha*x - alpha*alpha*x*x */
314  *this = Interval<Treal>::intersect(*this,oneInt);
315  // assert(upperBound <= 1.0);
316  upperBound = 2 * alpha * upperBound - alpha * alpha * upperBound * upperBound;
317  lowerBound = 2 * alpha * lowerBound - alpha * alpha * lowerBound * lowerBound;
318  if(nonEmptyIntervalInZeroToOne && upperBound < lowerBound) {
319  // Force interval to be non-empty.
320  Treal midPoint = (lowerBound + upperBound) / 2;
321  upperBound = lowerBound = midPoint;
322  }
323  }
324  else { /* ( alpha*x + (1-alpha) )^2 */
325  *this = Interval<Treal>::intersect(*this,zeroInt);
326  // assert(lowerBound >= 0.0);
327  upperBound = ( alpha * upperBound + (1-alpha) ) * ( alpha * upperBound + (1-alpha) );
328  lowerBound = ( alpha * lowerBound + (1-alpha) ) * ( alpha * lowerBound + (1-alpha) );
329  }
330  }
331 
332  template<typename Treal>
333  void Interval<Treal>::invPuriStep(int poly, Treal alpha) {
334  if (upperBound > 2.0 || lowerBound < -1.0)
335  throw Failure("Interval<Treal>::invPuriStep(int, real) : It is assumed here "
336  "that the interval I is within [-1.0, 2.0]");
337  Interval<Treal> oneInt(-1.0,1.0);
338  Interval<Treal> zeroInt(0.0,2.0);
339  if (poly) {
340  /* ( 1 - sqrt(1 - x) ) / alpha */
341  *this = Interval<Treal>::intersect(*this,oneInt);
342  // assert(upperBound <= 1.0);
343  upperBound = ( 1.0 - template_blas_sqrt(1.0 - upperBound) ) / alpha;
344  lowerBound = ( 1.0 - template_blas_sqrt(1.0 - lowerBound) ) / alpha;
345  }
346  else { /* ( sqrt(x) - (1-alpha) ) / alpha */
347  *this = Interval<Treal>::intersect(*this,zeroInt);
348  // assert(lowerBound >= 0.0);
349  upperBound = ( template_blas_sqrt(upperBound) - (1-alpha) ) / alpha;
350  lowerBound = ( template_blas_sqrt(lowerBound) - (1-alpha) ) / alpha;
351  }
352  }
353 
354 #endif
355 
356  template<typename Treal>
357  std::ostream& operator<<(std::ostream& s,
358  Interval<Treal> const & in) {
359  if (in.empty())
360  s<<"[empty]";
361  else
362  s<<"["<<in.low()<<", "<<in.upp()<<"]";
363  return s;
364  }
365 
366 } /* end namespace mat */
367 #endif
#define A
Interval< Treal > operator*(Treal const &value) const
Definition: Interval.h:154
bool overlap(Interval const &other) const
Definition: Interval.h:123
Treal upp() const
Definition: Interval.h:143
Treal length() const
Returns the length of the interval.
Definition: Interval.h:107
void decrease(Treal const value)
Definition: Interval.h:138
Interval(Treal low=1, Treal upp=-1)
Definition: Interval.h:46
Interval< Treal > operator+(Treal const &value) const
Definition: Interval.h:181
static Interval intersect(Interval const &A, Interval const &B)
Definition: Interval.h:51
Interval< Treal > operator/(Treal const &value) const
Definition: Interval.h:172
void increase(Treal const value)
Increases interval with value in both directions.
Definition: Interval.h:131
bool empty() const
Definition: Interval.h:49
Definition: allocate.cc:30
Treal midPoint() const
Definition: Interval.h:113
Definition: Interval.h:44
void invPuriStep(int poly)
Definition: Interval.h:269
void puriStep(int poly)
Definition: Interval.h:228
XYZpUV< TX, TY, TZ, TU, TV > operator+(XYZ< TX, TY, TZ > const &ABC, XY< TU, TV > const &DE)
Addition of two multiplication proxys XYZ and XY.
Definition: matrix_proxy.h:235
Interval< Treal > operator-(Treal const &value) const
Definition: Interval.h:178
Treal low() const
Definition: Interval.h:142
bool cover(Treal const value) const
Definition: Interval.h:117
Interval< Treal > sqrtInt(Interval< Treal > const &other)
Definition: Interval.h:217
Treal upperBound
Definition: Interval.h:195
std::ostream & operator<<(std::ostream &s, Interval< Treal > const &in)
Definition: Interval.h:357
void intersect(Interval const &other)
Definition: Interval.h:62
Interval< Treal > operator+(Interval< Treal > const &other) const
Definition: Interval.h:168
Definition: Failure.h:47
Treal lowerBound
Definition: Interval.h:194
Interval< Treal > operator-(Interval< Treal > const &other) const
Definition: Interval.h:162
XmY< TX, TY > operator-(TX const &AA, TY const &BB)
Substraction of two objects of type TX and TY.
Definition: matrix_proxy.h:275
#define B
void intersect_always_non_empty(Interval const &other)
Definition: Interval.h:78
Treal template_blas_sqrt(Treal x)