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