ergo
Interval.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 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>
42namespace 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
81 if (this->empty()) {
82 *this = other;
83 return;
84 }
85 if (other.empty())
86 return;
87
89 if ( !intersection.empty() ) {
90 *this = intersection;
91 return;
92 }
93 // Ok, the intersection is actually empty
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
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? */
165 return *this + (-1.0 * other);
166 // return Interval<Treal>(lowerBound - other.lowerBound,
167 // upperBound - other.upperBound);
168 }
169
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:
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);
225 }
226 }
227#endif
228
229 template<typename Treal>
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)
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)
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
Definition Failure.h:57
Definition Interval.h:46
void puriStep(int poly)
Definition Interval.h:230
void puriStep(int poly, Treal alpha)
Definition Interval.h:294
Treal low() const
Definition Interval.h:144
bool overlap(Interval const &other) const
Definition Interval.h:125
Interval< Treal > operator-(Treal const &value) const
Definition Interval.h:180
void increase(Treal const value)
Increases interval with value in both directions.
Definition Interval.h:133
Interval< Treal > operator/(Treal const &value) const
Definition Interval.h:174
void intersect_always_non_empty(Interval const &other)
Definition Interval.h:80
bool cover(Treal const value) const
Definition Interval.h:119
Treal midPoint() const
Definition Interval.h:115
Treal length() const
Returns the length of the interval.
Definition Interval.h:109
static Interval intersect(Interval const &A, Interval const &B)
Definition Interval.h:53
Treal upperBound
Definition Interval.h:197
Treal upp() const
Definition Interval.h:145
bool empty() const
Definition Interval.h:51
Interval< Treal > operator+(Treal const &value) const
Definition Interval.h:183
void invPuriStep(int poly)
Definition Interval.h:271
void decrease(Treal const value)
Definition Interval.h:140
void intersect(Interval const &other)
Definition Interval.h:64
Treal lowerBound
Definition Interval.h:196
Interval< Treal > operator*(Treal const &value) const
Definition Interval.h:156
void invPuriStep(int poly, Treal alpha)
Definition Interval.h:335
Interval< Treal > operator+(Interval< Treal > const &other) const
Definition Interval.h:170
Interval(Treal low=1, Treal upp=-1)
Definition Interval.h:48
Interval< Treal > operator-(Interval< Treal > const &other) const
Definition Interval.h:164
#define B
#define A
Definition allocate.cc:39
XmY< TX, TY > operator-(TX const &AA, TY const &BB)
Substraction of two objects of type TX and TY.
Definition matrix_proxy.h:277
Interval< Treal > sqrtInt(Interval< Treal > const &other)
Definition Interval.h:219
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
static Treal getMachineEpsilon()
Definition matInclude.h:147
std::ostream & operator<<(std::ostream &s, Interval< Treal > const &in)
Definition Interval.h:359
Treal template_blas_sqrt(Treal x)